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ABSTRACT 

The standard general relativistic model of a razor-thin accretion disk around a black hole, 
developed by Novikov & Thorne (NT) in 1973, assumes the shear stress vanishes at the radius 
of the innermost stable circular orbit (ISCO) and that, outside the ISCO, the shear stress is 
produced by an effective turbulent viscosity. However, astrophysical accretion disks are not 
razor-thin, it is uncertain whether the shear stress necessarily vanishes at the ISCO, and the 
magnetic field, which is thought to drive turbulence in disks, may contain large-scale struc- 
tures that do not behave like a simple local scalar viscosity. We describe three-dimensional 
general relativistic magnetohydrodynamic simulations of accretion disks around black holes 
with a range of spin parameters, and we use the simulations to assess the validity of the NT 
model. Our fiducial initial magnetic field consists of multiple (alternating polarity) poloidal 
field loops whose shape is roughly isotropic in the disk in order to match the isotropic turbu- 
lence expected in the poloidal plane. For a thin disk with an aspect ratio \hlr\ ~ 0.07 around 
a non-spinning black hole, we find a decrease in the accreted specific angular momentum of 
2.9% relative to the NT model and an excess luminosity from inside the ISCO of 3.5%. The 
deviations in the case of spinning black holes are also of the same order. In addition, the de- 
viations decrease with decreasing \h/r\. We therefore conclude that magnetized thin accretion 
disks in x-ray binaries in the thermal/high-soft spectral state ought to be well-described by the 
NT model, especially at luminosities below 30% of Eddington where we expect a very small 
disk thickness \h/r\ < 0.05 . We use our results to determine the spin equilibrium of black hole 
accretion disks with a range of thicknesses and to determine how electromagnetic stresses 
within the ISCO depend upon black hole spin and disk thickness. We find that the electro- 
magnetic stress and the luminosity inside the ISCO depend on the assumed initial magnetic 
field geometry. We consider a second geometry with field lines following density contours, 
which for thin disks leads to highly radially-elongated magnetic field lines. This gives roughly 
twice larger deviations from NT for both the accreted specific angular momentum and the lu- 
minosity inside the ISCO. Lastly, we find that the disk's corona (including any wind or jet) 
introduces deviations from NT in the specific angular momentum that are comparable to those 
contributed by the disk component, while the excess luminosity of bound gas from within the 
ISCO is dominated by only the disk component. Based on these indications, we suggest that 
differences in results between our work and other similar work are due to differences in the 
assumed initial magnetic field geometry as well as the inclusion of disk gas versus all the gas 
when comparing the specific angular momentum from the simulations with the NT model. 

Key words: accretion, accretion discs, black hole physics, hydrodynamics, (magnetohydro- 
dynamics) MHD, methods; numerical, gravitation 
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1 INTRODUCTION 

Accreting black holes (BHs) are among the most powerful astro- 
physical objects in the Universe. Although they have been the tar- 
get of intense research for a long time, many aspects of black hole 
accretion theory remain uncertain to t h is day . Pioneering work by 
iBarde en (l91U); ' Shakura & SunvaevI ( Il973l) : iNovikov & Thomel 
( ll973l) : iPage & Thomd h974h and"others indicated that black hole 
accretion through a razor-thin disk can be highly efficient, with up 
to 42% of the accreted rest-mass-energy being converted into ra- 
diation. These authors postulated the existence of a turbulent vis- 
cosity in the disk, parame terized via the famous ^-prescription 
dShakura & Sunvae\j[l973h . This viscosity causes outward trans- 
port of angular momentum; in the process, it dissipates energy and 
produces the radiation. The authors also assumed that, within the 
inner-most stable circular orbit (ISCO) of the black hole, the vis- 
cous torque vanishes and material plunges into the black hole with 
constant energy and angular momentum flux per unit rest-mass 
flux. This is the so-called "zero-torque" boundary condition. 

Modern viscous hydrodynamical calculations of disks with 
arbitrary thicknesses suggest that the zero-torque condition is 
a good approximation when the height (h) t o radius (r ) ratio 
of th e accreting gas is small: \h/r\ < 0.1 (Paczvhski 
Afshordi & Paczvnski 2003; Shafee et al. 2008b; Sadowski 



Abramowicz et all201(]|) . Radiatively efficient disks in active galac- 



tic nuclei (AGN) and x-ray binaries are expected to have disk thick- 
ness \h/r\ < 0.1 whenever the luminosi ty is limited to less tha n 
about 30% of the Eddington luminosity iMcClintock et alj|200^ . 
The above hydrodynamical studies thus suggest that systems in this 
limit should be described well by the standard relativ i stic th in disk 
theory as originally developed bv lNovikov & Thomel l fl973h . here- 
after NT. 

In parallel with the above work, it has for long been rec- 
ognized that the magnetic field could be a complicating factor 
that may s ignificantly m odify accretion dynamics near and inside 
the ISCO l lThomelll974]) . This issue has become increasingly im- 
portant with the realization that angular momentum transport in 
disks is entirely due to tu rbulence generated via the mag netoro- 
tational instability (MRI) jBalbus & Hawlevlll99ll. Il998h . How- 
ever, the magnetic field does not necessarily behave like a local 
viscous hydrodynamical stress. Near t he black hole, the m agnetic 
field may have large-scale structures I MacDonal j 1984). which 
can induce stresses across the ISCO jKrolikll 19991 ; lGammie|[l999l ; 
lAgol & Kroli3l200(ll) leading to changes in, e.g., the density, ve- 
locity, and amount of dissipation and emission. Unlike turbulence, 
the magnetic fie ld can transport angular momentum without dis- 
sipation (e.g. lLil l2002h . or it can dissip ate in current sheets w ith- 
out transporting angular momentum. In lAgol & KrolikI | |2000|) . the 
additional electromagnetic stresses are treated simply as a freely 
tunable model parameter on top of an otherwise hydrodynamical 
model. A more complete magnetohydrodynamic al (MHD ) mode l 
of a magnetized thin disk has been developed bv lGammiel l ll999l) . 
In this model, the controlling free parameter is the specific mag- 
netic flux, i.e., magnetic flux per unit rest-mass flux. Larger values 
of this parameter lead to larger deviations from NT due to electro- 
magnetic stresses, but the exact value of the parameter for a given 
accretion disk is unknown. For instance, it is entirely possible that 
electromagnetic stresses become negligible in the limit when the 
disk thickness \h/r\ — » 0. The value of the specific magnetic flux 
is determined by the nonlinear turbulent saturation of the magnetic 
field, so accretion disk simulations are the best way to establish its 
magnitude. 



The coupling via the magnetic field between a spinning 
black hole and an external disk, or between the hole and the 
corona, wind and jet (hereafter, corona-wind-jet), might also 
play an important role in modifying the accretion flow near 
the black hole. The wind or jet (hereafter, wind-jet) can trans- 
port angular momentum and energy awa y from the accretion 
disk an d black hole (iBlandfordI Il976l ; iBlandford & Znaiekl 
19771; iMcKinnev & Gammid |2004 iMcKinnevI 12006^ ; 



McKinnev & Naravaj l2007bl ; iKomissarov & McKinnevI |2007|) . 
The wind-jet power depen ds upon factors such as the black hole 
pin tMcKinnev 2005 ; Hawley & KrolikI 20061; Komissaroy et al .1 
2007I) . disk thickness jMeied l200ll ; iTchekhovskov et all l200a 



2009al lbl bold) , an d the strength and large-sc a le behavior of 



the m a gnetic field i McKinnev & Gammid 1 20041 ; iBeckwith et al.l 
l2008al ; iMcKinnev & Blandfor j |2009|) . and these can afi'ect the 
angular momentum transport through an accretion disk. In this 
context, we note that understanding how such factors affect disk 
structure may be key in interpreting the distinct states and complex 
behaviors observed for black hole X -ray binaries (Fender et al] 
|2004 iRemillard & McClintockl |2006|). These f actors also affect 
the black hole spin history jGammie et ai]|2004l) . and so must be 
taken into account when considering the eff^e ct of accretion on the 
cosmological evolution of black hole spin (Hu ghes & BlandfordI 
[2003; Gammie et al. 2004; Berti & Volonteri 20(M- 

Global simulations of accretion disks using ge neral relativis- 
tic magnetohydrodynamical (GRMHD) codes (e.g. iGammie et all 
l2003l ; lDe Villiers et al. Illooi currently provide the most complete 
understanding of how turbulent magnetized accretion flows around 
black holes work. Most simulations have studied thick (\h/r\ > 

0. 15) disks without radiative cooling. Such global simulations of 
the inner accretion flow have shown that fluid crosses the ISCO 
without any clear evidence that the torque vanishes at the ISCO, 

1. e., there is no apparent "stress edge" dMcKinnev & Gammid 
|2004 iKrolik et aljHoOSl ; IBeckwith et alj l2008hir Similar results 
were previ ously found with a pseu do-Newtonian potential for the 
black hole dKrolik & Hawlevll2002l) . In these studies, a plot of the 
radial profile of the normalized stress within the ISCO appears 
to indicate a significant deviation from th e NT thin disk theory 
dKrolik et al] |2005| ; ISeckwith et all l2008bl) . and it was thus ex- 
pected that much thinner disks might also deviate significantly from 
NT. A complicating factor in drawing firm conclusions from such 
studies is that the assumed initial global magnetic field geometry 
and strength can significantly change the magnitude of electromag- 
netic stres ses and associated angular momentum transport in side 
the ISCO dMcKinnev & Gammiell2004l ;l Beckwith et al ]|2008d) . 

The implications of the above studies for truly thin disks 
(\h/r\ < 0.1) remain uncertain. Thin disks are difficult to re- 
solve numerically, and simulations have been attempted only re- 
cently. Simulations of thin disks using a pseudo-Newtonian poten- 
tial for t he black hole reveal good a greement with standard thin disk 
theory dRevnolds & FabianlHoO^ . The first simulation of a thin 
(\h/r\ « 0.05) disk using a full GRMHD model was bv lShafee et all 
d2008d) . hereafter S08. They considered a non-spinning (a/M = 0) 
black hole and an initial field geometry consisting of multiple 
opposite-polarity bundles of poloidal loops within the disk. They 
found that, although the stress profile appears to indicate signifi- 
cant torques inside the ISCO, the actual angular momentum flux 
per unit rest-mass flux through the disk component deviates from 
the NT prediction by only 2%, corresponding to an estimated de- 
viation in the lum inosity of on l y abo ut 4%. The study by S08 was 
complemented bv lNobleetai] ( l2009h . hereafter N09, who consid- 
ered a thin (\h/r\ x 0.1) disk around an a/M = 0.9 black hole 
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and an initial field geometry consisting of a single highly-elongated 
poloidal loop bundle whose field lines follow the density contours 
of the thin disk. Th ey found 6% more l uminosity than predicted by 
NT. More recently, iNoble et alj feoiol) . hereafter NlO, considered 
a thin {\h/r\ x 0.07) disk around a non-spinning (a/M = 0) black 
hole and reported up to 10% deviations from NT in the specific 
angular momentum accreted through the accretion flow. 

In this paper, we extend the work of SOS by considering a 
range of black hole spins, disk thicknesses, field geometries, box 
sizes, numerical resolutions, etc. Our primary result is that we con- 
firm SOS, viz., geometrically thin disks are well-described by the 
NT model. We show that there are important differences between 
the dynamics of the gas in the disk and in the corona-wind-jet. In 
addition, we find that the torque and luminosity within the ISCO 
can be significantly affected by the geometry and strength of the 
initial magnetic field, a result that should be considered when com- 
paring simulation results to thin disk theory. In this context, we dis- 
cuss likely reasons for the apparently different conclusions reached 
by N09 and NlO. 

The equations we solve are given in ^ diagnostics are de- 
scribed in ^ and our numerical setup is described in ^ Results 
for our fiducial thin disk model for a non-rotating black hole are 
given in ^ and we summarize convergence studies in ^ Results 
for a variety of black hole spins and disk thicknesses are presented 
in ^and for thin disks with different magnetic field geometries 
and strengths in ^ We compare our results with previous studies 
in ^ discuss the implications of our results in 3101 and conclude 
with a summary of the salient points in 3111 



2 GOVERNING EQUATIONS 

The system of interest to us is a magnetized accretion disk 
around a rotating black hole. We write the black hole Kerr metric 
in Kerr-Schild (KS, horizon-pen etrating) coordinates jpont et all 
ll99Sl : |Papadopoulos & Fontll99Sh , which can be mapped to Boyer- 
Lindquist (BL) coordinates o r an orthonormal basis in any frame 
jMcKinnev & Gammiell200i) . We work with Heaviside-Lorentz 
units, set the speed of light and gravitational constant to unity 
(c = G = 1), and let M be the black hole mass. We solve the 
general relativistic magnetohydrody namical (GRMHD) e quations 
of motion for rotating black holes dGammie etallllOOal) with an 
additional cooling term designed to keep the simulated accretion 
disk at a desired thickness. 
Mass conservation gives 



V^o"") = 0, 



(1) 



where po is the rest-mass density, corresponding to the mass density 
in the fluid frame, and if is the contravariant 4-velocity. Note that 
we write the orthonormal 3-velocity as Vj (the covariant 3-velocity 
is never used below). Energy-momentum conservation gives 



(2) 



where the stress energy tensor T, includes both matter and electro- 
magnetic terms, 

T'- = (po + Ug + pg + b^)u"u, + (pg + b^lTjS", - If by, (3) 

where Ug is the internal energy density and pg = (T - \)iig is the 
ideal gas pressure with F = 4/3 Q The contravariant fluid-frame 



' Models wit h F = 5/3 show some minor differences compared to mode ls 
with F = 4/3 iMcKinnev & Gammiel2004l;lMignone & McKinne^^l2007l) . 



magnetic 4-field is given by N', and is related to the lab-frame 3- 
field via = B''h'^,/u' where /j^J = u''uy + 6y is a projection tensor, 
and S'!, is the Kronecker delta function. We write the orthonormal 3- 
field as S, (the covariant 3-field is never used below). The magnetic 
energy density (wj) and magnetic pressure (pj) are then given by 
"i = Pb = l^bf,/! = b- 12. Note that the angular velocity of the gas 
is Q = ii'^ lu' . Equation ^ has a source term 

which is a radiation 4-force corresponding to a simple isotropic 
comoving cooling term given by dU IdT. We ignore radiative trans- 
port effects such as heat currents, viscous stresses, or other effects 
that would enter as additional momentum sources in the comoving 
frame. In order to keep the accretion disk thin, we employ the same 
ad hoc cooling function as in SOS: 



— = -Ug 5[6)], 

dr r,ooi 



(5) 



where t is the fluid proper time, K = p^/p^ is the entropy con- 
stant, A",. = 0.00069 is set to be the same entropy constant as 
the torus atmosphere and is the entropy constant we cool the 
disk towards, and Kq > is the entropy constant of the ini- 
tial toru|^ The gas cooling time is set to t^^oI = 2;r/f2K, where 
Q.K = {\IM)l[(alM) + {RIMf'^] is the Keplerian angular frequency 
and R = r sin f is the cylindrical radius (We consider variations in 
the cooling timescale in section \5J\ ). We use a shaping function 
given by the quantity S [6] = exp[-(e - Jr/2)^/(2(flnocooi)^)], where 
we set Snocooi - (0.1,0.3,0.45,0.45} for our sequence of models 
with target thickness of \h/r\ = (0.07, 0.1, 0.2, 0.3), although we 
note that the thickest model with target \h/r\ = 0.3 has no cool- 
ing turned on. The shaping function 5 [6] is used to avoid cooling 
in the low density funnel-jet region where the thermodynamics is 
not accurately evolved and where the gas is mostly unbound (see 
Figure [T3] in section [5^ . In addition, we set the cooling function 
dU/dr = if 1) the timestep, dt, obeys dt > TcooI, which en- 
sures that the cooling does not create negative entropy gas ; or 2) 
log(K/Kc) < 0, which ensures the gas is only cooled, never heated. 
Photon capture by the black hole is not included, so the luminos- 
ity based upon this cooling function is an upper limit for radiation 
from the disk. The above cooling function drives the specific en- 
tropy of the gas toward the reference specific entropy Kc- Since 
specific entropy always increases due to dissipation, this cooling 
function explicitly tracks dissipation. Hence, the luminosity gen- 
erated from the cooling function should not be considered as the 
true luminosity, but instead should be considered as representing 
the emission rate in the limit that all dissipated energy is lost as 
radiation. Any other arbitrary cooling function that does not track 
dissipation would require full radiative transfer to obtain the true 
luminosity. 

^ We intended to have a constant K everywhere att = 0, but a normalization 
issue led to < Kg. Because of this condition, the disk cools toward a 
slightly thinner equilibrium at the start of the simulation, after which the 
cooling proceeds as originally desired by cooling towards the fiducial value 
K = Kc- Our models with |/i/r| x 0.07 are least affected by this issue. Also, 
since we do not make use of the cooling-based luminosity near ^ = 0, this 
issue does not affect any of our results. We confirmed that this change leads 
to no significant issues for either the magnitude or scaling of quantities with 
thickness by repeating some simulations with the intended = Kq. The 
otherwise similar simulations have thicker disks as expected (very minor 
change for our thin disk model as expected), and we find consistent results 
for a given measured thickness in the saturated state. 
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Magnetic flux conservation is given by tiie induction equation 

a,{ ^f^B-) = -dj[^r-g{B-v' - B'v-)], (6) 

wiiere v' = u' jii' , and g = Det(g,,y) is the determinant of the met- 
ric. In steady-state, the cooling is balanced by heating from shocks, 
grid-scale reconnection, and grid-scale viscosity. No explicit resis- 
tivity or viscosity is included. 



3 DIAGNOSTICS 

In this section, we describe several important diagnostics that we 
have found useful in this study. First, since we regulate the disk 
height via an ad hoc cooling function, we check the scale height 
of the simulated disk as a function of time and radius to ensure 
that our cooling function operates properly. Second, the equations 
we solve consist of seven time-dependent ideal MHD equations, 
corresponding to four relevant conserved quantitiefl Using these 
quantities we construct three dimensionless flux ratios correspond- 
ing to the accreted specific energy, specific angular momentum, and 
specific magnetic flux. Third, we check what the duration of the 
simulations should be in order to reach a quasi-steady state ("in- 
flow equilibrium") at any given radius. Finally, we describe how 
we compute the luminosity. 

When the specific fluxes are computed as a spatial or temporal 
average/integral, we always take the ratio of averages/integrals of 
fluxes (i.e. j d.xFi / j dxF2) rather than the average/integral of the 
ratio of fluxes (i.e. j dx{F[/ F2)). The former is more appropriate 
for capturing the mean behavior, while the latter can be more ap- 
propriate when investigating fluxes with significant phase shifted 
correlations between each other. As relevant for this study, the ac- 
cretion disk has significant vertical stratification and the local value 
of the ratio of fluxes can vary considerably without any affect on 
the bulk accretion flow. Similarly, potentially one flux can (e.g.) 
nearly vanish over short periods, while the other flux does not, 
which leads to unbounded values for the ratio of fluxes. However, 
the time-averaged behavior of the flow is not greatly affected by 
such short periods of time. These cases demonstrate why the ra- 
tio of averages/integrals is always performed for both spatial and 
temporal averages/integrals. 

When comparing the flux ratios or luminosities from the sim- 
ulations against the NT model, we evaluate the percent relative dif- 
ference D[f] between a quantity / and its NT value as follows: 

../■-./■[NT] 



D[f] = 100- 

/[NT] 

For a density-weighted time-averaged value of /, we compute 
/ / j fPo(r,e,(l>)dAg^dt 



</>P0 = 



/ / j'poir,9,cf>)dAg^dt 



(7) 



(8) 



where dAg^ = -\/^ddd(p is an area element in the 6-0 plane, and 
the integral over dt is a time average over the duration of interest, 
which corresponds to the period when the disk is in steady-state. 
For a surface-averaged value of /, we compute 



jjdAs^ 



(9) 



^ The energy-momentum of the fluid is not strictly conserved because of 
radiative cooling; however, the fluid component of the energy-momentum 
equations still proves to be useful. Only energy conservation of the fluid is 
strongly affected for our types of models. 
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3.1 Disk Thickness Measurement 

We define the dimensionless disk thickness per unit radius, \h/r\, 
as the density-weighted mean angular deviation of the gas from the 
midplane. 



(10) 



(This quantity was called AS^bs in SOS.) Notice that we assume the 
accretion disk plane is on the equator (i.e. we assume {0)pu = n/2). 
As defined above, \h/r\ is a function of r. When we wish to char- 
acterize the disk by means of a single estimate of its thickness, 
we use the value of \h/r\ at r = 2risco> where nsco is the ISCO 
radius (nsco = 6M for a non-spinning BH and nsc o = M for a 
maximally-spinning BH; IShapiro & Teukolskvll 19831) . As we show 
in ^5.41 this choice is quite reasonable. An alternative thickness 
measure, given by the root-mean-square thickness (/!/'')rms. allows 
us to estimate how accurate we can be about our definition of thick- 
ness. This quantity is defined by 



The range of 6 for the disk thickness integrals in the above equa- 
tions is from to n. 



3.2 Fluxes of Mass, Energy, and Angular Momentum 

The mass, energy and angular momentum conservation equations 
give the following fluxes. 



MO; t) 

E{r,t) 
M{r, t) 

J(r,l) 



nPQll'' 

Mir, t) ' 



dA„. 



(12) 



(13) 



(14) 



M(r, t) M(r, t) 

The above relations define the rest-mass accretion rate (sometimes 
just referred to as the mass accretion rate), M\ the accreted energy 
flux per unit rest-mass flux, or specific energy, e; and the accreted 
angular momentum flux per unit rest-mass flux, or specific angular 
momentum, j. Positive values of these quantities correspond to an 
inward flux through the black hole horizon. 

The black hole spin evolves due to the accretion of mass, en- 
ergy, and angular momentum, which can be described by the di- 
mensionless spin-up parameter s, 

_ dia/M) M 



dt M 



J 



M 



(15) 



where the ang ular integrals used to compute j and e include all 9 
and (p angles taammieetai]|2004l) . For s = the black hole is 
in so-called "spin equilibrium," corresponding to when the dimen- 
sionless black hole spin, a/M, does not change in time. 

The "nominal" efficiency, corresponding to the total loss of 
specific energy from the fluid, is obtained by removing the rest- 
mass term from the accreted specific energy: 



e = 1 - e. 



(16) 



The time-averaged value of e at the horizon (r = th) gives the total 
nominal efficiency: <e(rH)), which is an upper bound on the total 
photon radiative efficiency. 

The range of 9 over which the flux density integrals in the 
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above equations are computed depends on the situation. In SOS, 
we limited the 9 range to 69 = ±0.2 corresponding to 2-3 density 
scale heights in order to focus on the disk and to avoid including 
the disk wind or black hole jet. In this paper, we are interested in 
studying how the contributions to the fluxes vary as a function of 
height above the equatorial plane. Our expectation is that the disk 
and corona-wind-jet contribute differently to these fluxes. Thus, we 
consider different ranges of 9 in the integrals, e.g., from (;r/2) - 
2\h/r\ to (;r/2) + 2\h/r\, (tt/I) - 4\h/r\ to (;r/2) -I- 4\h/r\, or to n. 
The first and third are most often used in later sections. 



3.3 Splitting Angular Momentum Flux into Ingoing and 
Outgoing Components 

For a more detailed comparison of the simulation results with the 
NT model, we decompose the flux of angular momentum into an in- 
going ("in") term which is related to the advection of mass-energy 
into the black hole and an outgoing ("out") term which is related 
to the forces and stresses that transport angular momentum radially 
outward. These ingoing and outgoing components of the specific 
angular momentum are defined by 



Jm(r, t) 



((fio + iig + b-/2)u''){u^} 



{-Poll') 
7out('-,0 = J-Jm{r,t). 



(17) 
(18) 



By this definition, the "in" quantities correspond to inward trans- 
port of the comoving mass-energy density of the disk, ifiC'T^y = 
po + Ug + b- 12. Note that "in" quantities are products of the mean 
velocity fields {u' ) and (;(^) and not the combination (m' m^j); the lat- 
ter includes a contribution from correlated fluctuations in u!" and u^, 
which corresponds to the Reynolds stress. The residual of the total 
flux minus the "in" flux gives the outward, mechanical transport 
by Reynolds stresses and electromagnetic stresses. One could also 
consider a similar splitting for the specific energy. The above de- 
composition most closely matches our expectation that the inward 
flux should agree with the NT result as — > 0. Note, however, 
that our conclusions in this paper do not require any particular de- 
composition. This decomposition is different from SOS and NIO 
where the entire magnetic term (b^u' u^ - b'b^) is designated as the 
"out" term. Their choice overestimates the effect of electromag- 
netic stresses, since some magnetic energy is simply advected into 
the black hole. Also, the splitting used in SOS gives non-monotonic 
Jin vs. radius for some black hole spins, while the splitting we use 
gives monotonic values for all black hole spins. 



3.4 The Magnetic Flux 

The no-monopoles constraint implies that the total magnetic flux 
(O = §-dA) vanishes through any closed surface or any open sur- 
face penetrating a bounded flux bundle. The magnetic flux conser- 
vation equations require that magnetic flux can only be transported 
to the black hole or through the equatorial plane by advection. The 
absolute magnetic flux (J^, \B ■ dA\) has no such restrictions and can 
grow arbitrarily due to the MRI. However, the absolute magnetic 
flux can saturate when the electromagnetic field comes into force 
balance with the matter. We are interested in such a saturated state 
of the magnetic field within the accretion flow and threading the 
black hole. 

We consider the absolute magnetic flux penetrating a spherical 



surface and an equatorial surface given, respectively, by 

<b,ir,9j) = f f \B'\dAg,^, 
Je J4, 

<fo(r,9,t) = I I {BVA,,.. 



17' 



(19) 



(20) 



Nominally, O,- has an integration range of 9' = to 0' = 9 when 
measured on the black hole horizon, while when computing quanti- 
ties around the equatorial plane ff has the range {0) ± 9. One useful 
normalization of the magnetic fluxes is by the total flux through one 
hemisphere of the black hole plus through the equator 

(I>,o,(r, f) = 0,(r' = th, e' = . . . 7r/2, f) + (b^ir, ff = ;r/2, f), (21) 

which gives the normalized absolute radial magnetic flux 

<br{r, 9, t) 



*, ('•, e, t) ; 



(22) 



rt(r = .Rou,, ? = 0) ' 

where i?o„, is the outer radius of the computational box. The nor- 
malized absolute magnetic flux measures the absolute magnetic 
flux on the black hole horizon or radially through the equatorial 
disk per unit absolute flux that is initially available. 

The lGammid ( 1 19991) model of a magnetized thin accretion 
flow suggests another useful normalization of the magnetic flux is 
by the absolute mass accretion rate 



Mair, t) 



= { { Pol"'' 



(23) 



which gives the normalized specific absolute magnetic fluxes 
<t>r(r, t) 



E(r,t) 



Y(r,0 



Mcir, t) 
H V2 



(24) 





i 


Ma(r = rn,l) 


M 




SAh 



(25) 



where SA = (l/r^) j^^dAg^ is the local solid angle, SAh = 
SA(r = th) is the local solid angle on the horizon, H(r, t) is the ra- 
dial magnetic flux per unit rest-mass flux (usually specific magnetic 
flux), and T(r, t)cr'^^ jG is a particular dimensionless normalization 
of the specific mag netic flux that ap pears in the MHD accretion 
model developed bv lGamrni3 ( Il999l) . Since the units used for the 
magnetic field are arbitrary, any constant factor can be applied to 
H and one would still identify the quantity as the specific magnetic 
flux. So to simplify the discussion we henceforth call T the spe- 
cific magnetic flux. To obtain Equation ( I25t . all involved integrals 
should have a common 9 range around the equator. These quanti- 
ties all have absolute magnitudes because a sign change does not 
change the physical effect. The quantities 7, e, e, E, and T are each 
conserved along poloidal field-flow lines for stationary ide al MHD 
solutions l lBekenstein&Oronlll97Sl : iTakahashi et aDll990h . 

Gammie's (1999) model of a magnetized accretion flow within 
the ISCO assumes: 1) a thin equatorial flow ; 2) a radial poloidal 
field geometry (i.e., \Bg\ <K \B,\) ; 3) a boundary condition at the 
ISCO corresponding to zero radial velocity ; and 4) no thermal con- 
tribution. The model reduces to the NT solution within the ISCO 
for Y — » 0, and deviations from NT's solution are typically small 
(less than 12% for 7 across all black hole spins; see Appendix [AJ 
for T < 1 . We have defined the quantity T in equation ( 124b with the 
V2 factor, the square root of the total mass accretion rate through 
the horizon per unit solid angle, and Heaviside-Lorentz units for B' 
so that the numerical value of T at the horizon is identical l y equa l 
to the numerical value of the free parameter in lOammi 3 ll999h . 
i.e., their Fg^ normalized by Fm = -1. As shown in that paper. 
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T directly controls deviations of the specific angular momentum 
and specific energy away from the non-magnetized thin disk the- 
ory values of the NT model. Even for disks of finite thickness, 
the parameter shows how electromagnetic stresses control devia- 
tions between the horizon and the ISCO. Note that the flow can 
devia te from NT at the ISCO sim ply due to finite thermal pres- 
sure jMcKinnev & Gamrni3l2004l) . In Appendix I A II Table El we 
list numerical values of j and e for Gammie's (1999) model, and 
show how these quantities deviate from NT for a given black hole 
spin and Y. 

We find T to be more useful as a measure of the importance of 
the magnetic field within the ISCO than our previous measurement 
in SOS of the a-viscosity parameter, 

a = , (26) 

Ps + Pb 

where T'^'' = e*,,e',T'"' is the orthonormal stress-energy tensor 
components in the comoving frame, and e^^, is the contravariant 
tetrad system in the local fluid-frame. This is related to the normal- 
ized stress by 



M 



(27) 



where dA'g^ = e^^e'^ydff^dif)^' is the comoving area element, dL'^ = 
e^yd(fi^' evaluated at 9 = n/2 is the comoving ip length element, 
ff" = {0,0, 1,0), and 0'' = {0,0,0, 1). This form for W is a simple 
generalization of Eq. 5.6.1b in NT73, and note that the NT solution 
for W/M is given by Eq. 5.6.14a in NT73. In SOS, W was inte- 
grated over fluid satisfying —u,{po + Ug + Pg + b^)/po < 1 (i.e., 
only approximately gravitationally bound fluid and no wind-jet). 
We use the same definition of bound in this paper. As shown in 
SOS, a plot of the radial profile of W/M or a within the ISCO does 
not necessarily quantify how much the magnetic field affects the 
accretion flow properties, since even apparently large values of this 
quantity within the ISCO do not cause a significant deviation from 
N T in the specific a ngular momentum accreted. On the other hand, 
the lGammiel ( ll999l) parameter Y does directly relate to the electro- 
magnetic stresses within the ISCO and is an ideal MHD invariant 
(so constant vs. radius) for a stationary flow. One expects that ap- 
propriately time-averaged simulation data could be framed in the 
context of this stationary model in order to measure the effects of 
electromagnetic stresses. 

3.5 Inflow Equilibrium 

When the accretion flow has achieved steady-state inside a given ra- 
dius, the quantity M(r, t) will (apart from normal fluctuations due to 
turbulence) be independent of time, and if it is integrated over all 6 
angles will be constant within the given radiu^. The energy and an- 
gular momentum fluxes have a non-conservative contribution due 
to the cooling function and therefore are not strictly constant. How- 
ever, the cooling is generally a minor contribution (especially in the 
case of the angular momentum flux), and so we may still measure 
the non-radiative terms to verify inflow equilibrium. 

The radius out to which inflow equilibrium can be achieved 
in a given time can be estimated by calculating the mean radial 

^ If we integrate over a restricted range of 9, then there could be addi- 
tional mass flow through the boundaries in the 6 direction and M(r, t) will 
no longer be independent of r, though it would still be independent of t. 



velocity i', and then deriving from it a viscous timescale -r/v,-. 
From standard accretion disk theory and using the definition of a 
given in Eq. l |27t , the mean radial velocity is given by 



(2S) 



where Wk ~ (r/M) ''^ is the Keplerian speed at radius r and a is the 
standard viscosity parameter given by equation l l26b jFrank et"aD 
[1992). Although the viscous timescale is the nominal time needed 
to achieve steady-state, in practice it takes several viscous times be- 
fo re the flow real ly settles down, e.g., see the calculations reported 
in IShapirdJ2O10l) ." In the present paper, we assume that inflow equi- 
librium is reached after two viscous times, and hence we estimate 
the inflow equilibrium time, /je, to be 



\3/2 



Vr \mI \a\h/rp 



I 



3/2 



5000 — I 



(29) 



where, in the right-most relation, we have taken a typical value of 
Q- ~ 0. 1 for the gas in the disk proper (i.e., outside the ISCO) and we 
have set \h/r\ a 0.064, as appropriate for our thinnest disk models. 

A simulation must run until t ~ fie before we can expect in- 
flow equilibrium at radius r. According to the above Newtonian 
estimate, a thin disk simulation with \hlr\ ~ 0.064 that has been 
run for a time of 30000M will achieve steady-state out to a ra- 
dius of only ~ 3M. However, this estimate is inaccurate since it 
does not allow for the boundary condition on the flow at the ISCO. 
Both the boundary condition as well as the effects of GR are in- 
cluded in the formula for the radial velocity given in Eq. 5.9.8 of 
NT, which we present for completeness in Appendix |B] That more 
accurate result, which is what we use for all our plots and numer- 
ical estimates, shows that our thin disk simulations should reach 
inflow equilibrium within r/M = 9, 7, 5.5, 5, respectively, for 
a/M = 0, 0.7, 0.9, 0.9S. These estimates are roughly consistent 
with the radii out to which we have a constant j vs. radius in the 
simulations discussed in ^ 



3.6 Luminosity Measurement 

We measure the radiative luminosity of the accreting gas directly 
from the cooling function dU/dr. At a given radius, r, in the steady 
region of the flow, the luminosity per unit rest-mass accretion rate 
interior to that radius is given by 



L{< r) 

M{r, t) M{r, l)(lf 



(tf - ti) J,=,, Jr'=n,^ Jfl=0 J0 \ dT I 



(30) 

where dV„Jn^ = yp^dtdr' d9d(fi and the 4D integral goes from the 
initial time /, to the final time tf over which the simulation results 
are time-averaged, from the radius th of the horizon to the radius r 
of interest, and usually over all 9 and (/>. We find it useful to compare 
the simulations with thin disk theory by computing the ratio of the 
luminosity emitted inside the ISCO (per unit rest-mass accretion 
rate) to the total radiative efficiency of the NT model: 



^in — 



L{< nsco) 



(31) 



Me[NT] 

This ratio measures the excess radiative luminosity from inside the 
ISCO in the simulation relative to the total luminosity in the NT 
model (which predicts zero luminosity here). We also consider the 
excess luminosity over the entire inflow equilibrium region 

_ L(r<req)-L(r<re,)[NT] 



Me[NT] 



(32) 
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which corresponds to the luminosity (per unit mass accretion rate) 
inside the inflow equilibrium region (i.e. r < r^^^, where is the 
radius out to which inflow equilibrium has been established) sub- 
tracted by the NT luminosity all divided by the total NT efficiency. 
Large percent values of Li„ or L^q would indicate large percent de- 
viations from NT. 



4 PHYSICAL MODELS AND NUMERICAL METHODS 

This section describes our physical models and numerical methods. 
Table [T] provides a list of all our simulations and shows the phys- 
ical and numerical parameters that we vary. Our primary models 
are labeled by names of the form AxHRy, where x is the value of 
the black hole spin parameter and y is approximately equal to the 
disk thickness \h/r\. For instance, our fiducial model A0HR07 has 
a non-spinning black hole (a/M = 0) and a geometrically thin disk 
with \h/r\ ~ 0.07. We discuss this particular model in detail in ^ 
Table [T] also shows the time span (from Tj/M to Tf/M) used to 
perform the time-averaging, and the last column shows the actual 
value of \h/r\ in the simulated model as measured during inflow 
equilibrium, e.g., \h/r\ = 0.064 for model A0HR07. 



4.1 Physical Models 

This study considers black hole accretion disk systems with a range 
of black hole spins: a/M = 0, 0.7, 0.9, 0.98, and a range of disk 
thicknesses: \h/r\ = 0.07, 0.13, 0.2, 0.3. The initial mass d is- 
tribution is an isentrop ic equilibrium torus l lChakrabartilll985j|bl 
iDe Villiers et alj2003l) . All models have an initial inner torus edge 
at Hn = 20M, while the torus pressure maximum for each model 
is located between R^.^^ = 35M and = 65M. We chose this 
relatively large radius for the initial torus because SOB found that 
placing the torus at smaller radii caused the results to be sensitive 
to the initial mass distribution. We initialize the solution so that 
po = 1 is the maximum rest-mass density. In SOS, we set q = 1.65 
(a oc r'' in non-relativistic limit) and K = 0.00034 with F = 4/3, 
while in this paper we adjust the initial angular momentum pro- 
file such that the initial torus has the target value of \h/r\ at the 
pressure maximum. For models with \h/r\ = 0.07, 0.13, 0.2, 0.3, 
we fix the specific entropy of the torus by setting, respectively, 
K = Ko = {0.00034, 0.0035, 0.009, 0.009) in the initial poly- 
tropic equation of state given by p = KqP^. The initial atmosphere 
surrounding the torus has the same polytropic equation of state with 
nearly the same entropy constant of = 0.00069, but with an ini- 
tial rest-mass density of po = 10"*(r/M)"''^, corresponding to a 
Bondi-like atmosphere. 

Recent GRMHD simulations of thick disks indicate that the 
results for the disk (but not the wind- jet) are roug hly indepen- 
dent of the initial fiel d geometry l lMcKinnev & Nara van 2007a, b; 
iBeckwith et al .l2008ah . However, a detailed study for thin disks has 
yet to be performed. 

We consider a range of magnetic field geometries described 
by the vector potential A,, which is related to the Faraday tensor by 



■ A^_y. As in SOS, we consider a general multiple-loop 



field geometry corresponding to N separate loop bundles stacked 
radially within the initial disk. The vector potential we use is given 
by 

, I log(r/S) \ 
oc Q- sin ^ I [1 + w(ranc - 0.5)] , (33) 



below for a discussion of perturbations.). All other A^ are initially 
zero. All our multi-loop and 1-loop simulations have S = 22M, 
and the values of Af^^t^iKlnr) are listed in Table [T] For multi-loop 
models, each additional field loop bundle has opposite polarity. We 
use Q = (j(j;/«g,n,ax - 0.2)(r/M)''''', and set 2 = if either r < S 
or 2 < 0, and Ug^^.^,^ is the maximum value of the internal energy 
density in the torus. By comparison, in SOS, we set 5 = Ll^n, Hn = 
20M, AfiM/(2nr) = 0.16, such that there were two loops centered 
at r = 2SM and 3SM. The intention of introducing multiple loop 
bundles is to keep the aspect ratio of the bundles roughly 1:1 in 
the poloidal plane, rather than loop(s) that are highly elongated in 
the radial direction. For each disk thickness, we tune Afj^t^/ilnr) in 
order to obtain initial poloidal loops that are roughly isotropic. 

As in SOS, the magnetic field strength is set such that the 
plasma /? parameter satisfies /?„i„es = Pg.max/Pft.max = 100, where 
Pg.max is the maximum thermal pressure and P6_,^ax is the maxi- 
mum magnetic pressure in the entire torus. Since the two max- 
ima never occur at the same location, p = Pg/pi, varies over 
a wide range of values within the disk. This approach is sim- 
ilar to how the magnetic field was normalized in other stud- 
ies taammie et alJbOoH; iMcKinnev & Gammi3l2004 iMcKinnevI 



McKinnev & Nara van 2007b ; K omissarov & McKinnevI 
It ensures that the magnetic field is weak throughout the 
disk. Care must be taken with how one normalizes any given initial 
magnetic field geometry. For example, f or the 1-loop field geome- 
try used bv lMcKinnev & Gammid ( |2004|) . if one initializes the field 
with a mean (volume-averaged) = 100, then the inner edge of the 
initial torus has /J ~ 1 and the initial disk is not weakly magnetized. 

For most models, the vector potential at all grid points was 
randomly perturbed by 2% (w in equation [33} and the internal en- 
ergy density at all grid points was randomly perturbed by 10% 0. 
If the simulation starts with perturbations of the vector potential, 
then we compute Oioi (used to obtain O,-) using the pre-perturbed 
magnetic flux in order to gauge how much flux is dissipated due 
to the perturbations. Perturbations should be large enough to ex- 
cite the non-axisymmetric MRI in order to avoid the axisymmetric 
channel solution, while they should not be so large as to induce sig- 
nificant dissipation of the magnetic energy due to grid-scale mag- 
netic dissipation just after the evolution begins. For some models, 
we studied different amplitudes for the initial perturbation in order 
to ensure that the amplitude does not significantly affect our re- 
sults. For a model with \h/r\ ~ 0.07, a/M = 0, and a single polarity 
field loop, one simulation was initialized with 2% vector potential 
perturbations and 10% internal energy perturbations, while another 
otherwise similar simulation was given no seed perturbations. Both 
become turbulent at about the same time t ~ 1500M. The mag- 
netic field energy at that time is negligibly different, and there is no 
evidence for significant differences in any quantities during inflow 
equilibrium. 



4.2 Numerical Methods 

We perform simulations using the GRMHD code HARM that 
is based upon a conservative shock-capturing Godunov scheme. 
One key feature of our code is that we us e horizon-penetrating 
Kerr-Schild coordinates for the Kerr metric i Gammie et al.| 2003 ; 
iMcKinnev & Gammiell2004l : lMcKinnevll2006j ; lNoble eTauboOe ; 



[ /lfleid/(2;Tr) j 

where ranc is a random number generator for the domain to 1 (see 



^ In SOS, we had a typo saying we perturbed the field by 50%, while it 
was actually perturbed the same as these models, i.e.: 2% vector potential 
perturbations and 10% internal energy perturbations. 
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Table 1. Simulation Parameters 
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iMignone & McKinnevI |2007| ; iTchekhovskov et alj |2007|) . which 
avoids any issues with the coordinate singularity in Boyer- 
Lindquist coordinates. Even with Kerr-Schild coordinates, one 
must ensure that the inner-radial boundary of the computational 
domain is outside the so-called inner horizon (at r/M = 1 - 
- (a/M)-) so that the equations remain hyperbolic, and one 
must ensure that there are plenty of grid cells spanning the region 
near the horizon in order to avoid numerical diffusion out of the 
horizon. 

Another key feature of our code is the use of a 3rd order ac- 
curate (4th order error) PPM scheme for the interpolation of primi- 
tive quantities (i.e. rest-mass density, 4- vel ocity relative to a Z AMO 
observer, and lab-frame 3-magnetic field) jMc Kinnev 2006b). The 
interpo lation is similar to that described in IColella & Woodward, 
( Il984l) . but we modified it to be consistent with interpolating 
through point values of primitives rather than average values. We 
do not use the PPM steepener, but we do use the PPM flattener that 
only activates in strong shocks (e.g. in the initial bow shock ofi" the 
torus surface, but rarely elsewhere). The PPM scheme attempts to 
fit a monotonic 3rd order polynomial directly through the grid face 
where the dissipative flux enters in the Godunov scheme. Only if 
the polynomial is non-monotonic does the interpolation reduce or- 
der and create discontinuities at the cell face, and so only then does 
it introduce dissipative fluxes. It therefore leads to extremely small 
dissipation compared to the original schemes used in HARM, such 
as the 1st order accurate (2nd order error) minmod or monotonized 
central (MC) limiter type schemes that always create discontinu- 
ities (and so dissipative fluxes) at the cell face regardless of the 
monotonicity for any primitive quantity that is not linear in space. 

Simulations of fully three-dimensional models of accreting 
black holes producing jets using our 3D GRMHD code show that 
this PPM scheme leads to an improvement in effective resolution 



by at least factors of roughly two per dimension as compared to 
the original HARM MC limiter s cheme for models with resolution 
256 X 128 X 32 (M cKinnev & Bl andford 2009). The PPM method 
is particularly well-suited for resolving turbulent flows since they 
rarely have strong discontinuities and have most of the turbulent 
power in long wavelength modes. Even moving discontinuities are 
much more accurately resolved by PPM than minmod or MC. For 
example, even without a steepener, a simple moving contact or 
moving magnetic rotational discontinuity is sharply resolved within 
about 4 cells using the PPM scheme as compared to being diffu- 
sively resolved within about 8-15 cells by the MC limiter scheme. 

A 2nd order Runge-Kutta method-of-lines scheme is used to 
step forward in time, and the timestep is set by using the fast mag- 
netosonic wavespeed with a Courant factor of 0.8. We found that a 
4th order Runge-Kutta scheme does not significantly improve ac- 
curacy, since most of the flow is far beyond the grid cells inside 
the horizon that determine the timestep. The standard HARM HLL 
scheme is used for the dissipative fluxes, and the standard HARM 
Toth scheme is used for the magnetic field evolution. 



4.3 Numerical Model Setup 

The code uses uniform internal coordinates (f, x*'', x*^', x*^') 
mapped to the physical coordinates (f, r, 9, (f>). The radial grid map- 
ping is 

r(x<") = i?o + exp(x"^), (34) 

which spans from Ri„ to i?out- The parameter Rq = 0.3M controls the 
resolution near the horizon. Absorbing (outflow, no inflow allowed) 
boundary conditions are used. The 6-grid mapping is 

9(x'->) = [y(2x'^' - 1) + (1 - y)(2x<^> - 1)' + l](7r/2), (35) 
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where x'^' ranges from to 1 (i.e. no cut-out at the poles) and Y is 
an adjustable parameter that can be used to concentrate grid zones 
toward the equator as Y is decreased from 1 to 0. Roughly half of 
the 6 resolution is concentrated in the disk region within ±2\h/r\ of 
the midplane. The HR07 and HR2 models listed in Table [T] have 
11 cells per \h/r\, while the HRl and HR3 models have 7 cells per 
\h/r\. The high resolution run, C6, has 22 cells per \h/r\, while the 
low resolution model, C5, has 5 cells per \h/r\. For y = 0.15 this 
grid gives ro ughly 6 times more angular res olution compared to the 
grid used in iMcKinnev & Gammid ( |2004) given by equation (8) 
with h = 0.3. Reflecting boundary conditions are used at the polar 
axes. 

The (/i-grid mapping is given by (pix^^^) = 27rx'", such that 

varies from to 1/8,1/4,3/8,1/2 for boxes with = 
n/4,n/2,3n/4,n, respectively. Periodic boundary conditions are 
used in the i/i-direction. In all cases, the spatial integrals are renor- 
malized to refer to the full 2;r range in (/>, even if our computa- 
tional box size is limited in the (/i-direction. We consider various 
A(f) in order to check whether this changes our results. Previous 
GRMHD simulations with the full A(f> = 27r extent suggest that 
A(fi = n/2 is sufficient since coh erent structures only ext end for 
about one radian (see Fig. 12 in ISchnittman et al] 1 20061) . How- 
ever, in other GRMHD studies with Aif> = 2n, the ni = 1 mode 
was found to be dominant, so this requires further consideration 
jMcKinnev & BlandfordI boOSl) . Note that SOS used A(/> = n/4, 
while both N09 and NIO used A(/> = n/2. 

The duration of our simulations with the thinnest disks varies 
from approximately 20000M to 30000M in order to reach inflow 
equilibrium and to minimize fluctuations in time-averaged quanti- 
ties. We ensure that each simulation runs for a couple of viscous 
times in order to reach inflow equilibrium over a reasonable range 
of radius. Note that the simulations cannot be run for a duration 
longer than ~ M^skit = 0)/M ~ 10^ M, corresponding to the 
time-scale for accreting a significant fraction of the initial torus. We 
are always well below this limit. 

Given finite computational resources, there is a competition 
between duration and resolution of a simulation. Our simulations 
run for relatively long durations, and we use a numerical resolution 
of N^y. NgX = 256 x 64 x 32 for all models (except those used 
for convergence testing). In S08 we found this resolution to be suf- 
ficient to obtain convergence compared to a similar 5 12 x 128 x 32 
model with A(p = n/4. In this paper, we explicitly confirm that 
our resolution is sufficient by convergence testing our results (see 
section [6]l. Near the equatorial plane at the ISCO, the grid aspect 
ratio in dr : rdd : rsinOdif) is 2:1:7, 1:1:4, 1:1:3, and 1:1:3, re- 
spectively, for our HR07, HRl, HR2, and HR3 models. The 2:1:7 
grid aspect ratio for the HR07 model was found to be sufficient in 
SOS. A grid aspect ratio of 1:1:1 would be preferable in order to 
ensure the dissipation is isotropic in Cartesian coordinates, since in 
Nature one would not expect highly anisotropic dissipation on the 
scale resolved by our grid cells. However, finite computational re- 
sources require a balance between a minimum required resolution, 
grid aspect ratio, and duration of the simulation. 

As described below, we ensure that the MRI is resolved in 
each simulation both as a function of space and as a function of 
time by measuring the number of grid cells per fastest growing MRI 
mode: 

_ _ /Imri , l<l/l^('-.^)l ,,,, 
Qum = — :— ~ 27T , iib) 

where Ag = \e''^,dxf \ is the comoving orthonormal 6-directed grid 



cell length, e*"^ is the contravariant tetrad system in the local fluid- 
frame, = -yjbgb^ l{b- + pa + Ug + pg) is the Alfven speed, b'' = 

e^i^V is the comoving orthonormal 0-directed 4-field, and \Q.{r, 0)\ 
is the temporally and azimuthally averaged absolute value of the 
orbital frequency. 

During the simulation, the rest-mass density and internal en- 
ergy densities can become quite low beyond the corona, but the 
code only remains accurate and stable for a finite value of b^/po, 
b-/ug, and Ug/po for any given resolution. We enforce b'/po < 10"*, 
b^/Ug < 10"*, and Ug/po < 10* by injecting a sufficient amount 
of mass or internal energy into a fixed zero angular momentum 
observer (ZAMO) frame with 4-velocity = j-a, 0, 0, 0), where 
a = 1/ y-g" is the lapse. In some simulations, we have to use 
stronger limits given by b-/pQ < 10, b^/iig < 10^, and iig/po g 10, 
in order to maintain stability and accuracy. Compared to our older 
method of injecting mass-energy into the comoving frame, the new 
method avoids run-away injection of energy-momentum in the low- 
density regions. We have confirmed that this procedure of injecting 
mass-energy does not contaminate our results for the accretion rates 
and other diagnostics. 



5 FIDUCIAL MODEL OF A THIN DISK AROUND A 
NON-ROTATING BLACK HOLE 

Our fiducial model, A0HR07, consists of a magnetized thin ac- 
cretion disk around a non-rotating (a/M = 0) black hole. This is 
similar to the model described in SOS; however, here we consider 
a larger suite of diagnostics, a resolution of 256 x 64 x 32, and a 
computational box with A(p = n/2. As mentioned in section|4T| the 
initial torus parameters are set so that the inner edge is at r = 20M, 
the pressure maximum is at r = 35M, and \h/r\ < 0.1 at the pressure 
maximum (see Figure[TJ. 

The initial torus is threaded with magnetic field in the multi- 
loop geometry as described in section |4T| For this model, we use 
four loops in order to ensure that the loops are roughly circular in 
the poloidal plane. Once the simulation begins, the MRI leads to 
MHD turbulence which causes angular momentum transport and 
drives the accretion flow to a quasi-steady state. 

The fiducial model is evolved for a total time of 27350M. We 
consider the period of steady-state to be from T, = 12500M to Tf = 
27350M and of duration AT = 14S50M. All the steady-state results 
described below are obtained by time-averaging quantities over this 
steady-state period, which corresponds to about 160 orbital periods 
at the ISCO, 26 orbits at the inner edge of the initial torus (r = 
20M), and 1 1 orbits at the pressure maximum of the initial torus 
{r = 35M). 



5.1 Initial and Evolved Disk Structure 

Figure[T]shows contour plots of various quantities in the initial so- 
lution projected on the (R, z) = (r sin 9, r cos 0)-plane. Notice the 
relatively small vertical extent of the torus. The disk has a thick- 
ness of \h /r\ ~ 0.06 - 0.09 over the radius range containing the bulk 
of the mass. The four magnetic loops are clearly delineated. The 
plot of 2mri indicates that the MRI is well-resolved within the two 
primary loops. The left-most and right-most loops are marginally 
under-resolved, so a slightly slower-growing MRI mode is expected 
to control the dynamics in this region. However, the two primary 
loops tend to dominate the overall evolution of the gas. 
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Figure |2] shows the time-averaged solution during the quasi- 
steady state period from T, = 12500M to Tf = 27350M. We refer 
to the disk during this period as being "evolved" or "saturated." The 
evolved disk is in steady-state up to r ~ 9M, as expected for the du- 
ration of our simulation. The rest-mass density is concentrated in 
the disk midplane within ±2\hlr\, while the magnetic energy den- 
sity is concentrated above the disk in a corona. The MRI is properly 
resolved with 2mri « 6 in the disk midplan^^. The gas in the mid- 
plane has plasma ~ 10 outside the ISCO and /3 ~ 1 near the black 
hole, indicating that the magnetic field has been amplified beyond 
the initial minimum of /? ~ 100. 

Figure [3] shows the time-averaged structure of the magnetic 
field during the quasi-steady state period. The field has a smooth 
split-monopole structure near and inside the ISCO. Beyond r ~ 
9M, however, the field becomes irregular, reversing direction more 
than once. At these radii, the simulation has not reached inflow 
equilibrium. 

5.2 Velocities and the Viscous Time-Scale 

Figure |4] shows the velocity structure in the evolved model. The 
snapshot indicates well-developed turbulence in the interior of the 
disk at radii beyond the ISCO (r > 6M), but laminar flow inside 
the ISCO and over most of the corona. The sudden transition from 
turbulent to laminar behavior at the ISCO, which is seen also in the 
magnetic field (Figure [3ji), is a clear sign that the flow dynamics 
are quite different in the two regions. Thus the ISCO clearly has 
an affect on the accreting gas. The time-averaged flow shows that 
turbulent fluctuations are smoothed out within r ~ 9M. Figure |5] 
shows the velocity stream lines using the line integral convolution 
method to illustrate vector fields. This figure again confirms that the 
accretion flow is turbulent at radii larger than nsco but it becomes 
laminar inside the ISCO, and it again shows that time-averaging 
smooths out turbulent fluctuations out to r ~ 9M. 

Figure|6]shows components of the time-averaged velocity that 
are angle-averaged over ±2\h/r\ around the midplane (thick dashed 
lines in Figure (S)). By limiting the range of the 6 integral, we fo- 
cus on the gas in the disk, leaving out the corona-wind-jet. Outside 
the ISCO, the radial velocity from the simulation agrees well with 
the analytical GR estimate (Eg. |B7| in Appendix iBt. By making 
this comparison, we found a\h/r\^ x 0.00033. For our disk thick- 
ness \h/r\ = 0.064, this corresponds to a a 0.08, which is slightly 
smaller than the nominal estimate a ~ 0. 1 we assumed in 33.51 As 
the gas approaches the ISCO, it accelerates rapidly in the radial di- 
rection and finally free-falls into the black hole. This region of the 
flow is not driven by viscosity and hence the dynamics here are not 
captured by the analytical formula. 

Figure [6] also shows the inflow equilibrium time fie, which we 
take to be twice the GR version of the viscous time: t^^ = -2r/v,.. 

ISano et alj i2004h find that having about 6 grid cells per wavelength of 
the fastest growing MRI mode during saturation leads to convergent behav- 
ior for the electromagnetic stresses, although their determination of 6 cells 
was based upon a 2nd order van Leer scheme that is significantly more dif- 
fusive than our PPM scheme. Also, the (time-averaged or single time value 
of) vertical field is already (at any random spatial position) paitially sheared 
by the axisymmetric MRI, and so may be less relevant than the (e.g.) max- 
imum vertical field per unit orbital time at any given point that is not yet 
sheared and so represents the vertical component one must resolve. These 
issues imply we may only need about 4 cells per wavelength of the fastest 
growing mode (as defined by using the time-averaged absolute vertical field 
strength). 




-6 ^^^^^^^^^^^^^^^^^^^^^"-1 

20 25 30 35 40 45 50 

r/M 



Figure 1. The initial state of the fiducial model (A0HR07) consists of 
weakly magnetized gas in a geometrically thin torus around a non- spinning 
(a/M = 0) black hole. Color maps have red as highest values and blue as 
lowest values. Panel (a): Linear color map of rest-mass density, with solid 
lines showing the thickness \h/r\ of the initial torus. Note that the black hole 
horizon is at r = 2M, fai' to the left of the plot, so the torus is clearly ge- 
ometrically thin. Near the pressure maximum \h/r\ < 0.1, and elsewhere 
\h/r\ is even smaller. Panel (b): Contour plot of overlaid on linear color 
map of rest-mass density shows that the initial field consists of four poloidal 
loops centered at r/M = 29, 34, 39, 45. The wiggles in are due to the 
initial perturbations. Panel (c): Linear color map of the plasma /? shows that 
the disk is weakly magnetized throughout the initial torus. Panel (d): Linear 
color map of the number of grid cells per fastest growing MRI wavelength, 
2mri, shows that the MRI is properly resolved for the primary two loops at 
the center of the disk. 

This is our estimate of the time it will take for the gas at a given 
radius to reach steady-state. We see that, in a time of ~ 27350M, 
the total duration of our simulation, the solution can be in steady- 
state only inside a radius of ~ 9M. Therefore, in the time-averaged 
results described below, we consider the results to be reliable only 
over this range of radius. 

5.3 Fluxes vs. Time 

Figure[7]shows various fluxes vs. time that should be roughly con- 
stant once inflow equilibrium has been reached. The figure shows 
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Figure 2. The evolved state of the fiducial model (A0HR07) consists of a 
weakly magnetized thin disk surrounded by a strongly magnetized corona. 
All plots show quantities that have been time-averaged over the period 
12500M to 27350M. Color maps have red as highest values and blue as 
lowest values. Panel (a): Linear color map of rest-mass density, with solid 
lines showing the disk thickness \h/r\. Note that the rest-mass density drops 
off rapidly inside the ISCO. Panel (b): Linear color map of shows that 
a strong magnetic field is present in the corona above the equatorial disk. 
Panel (c): Linear color map of plasma /J shows that the /? values are much 
lower than in the initial torus. This indicates that considerable field ampli- 
fication has occurred via the MRI. The gas near the equatorial plane has 
P ~ 10 far outside the ISCO and approaches /? ~ 1 near the black hole. 
Panel (d): Linear color map of the number of grid cells per fastest grow- 
ing MRI wavelength, 2mri. shows that the MRI is properly resolved within 
most of the accretion flow. Note that 2mri (determined by the vertical mag- 
netic field strength) is not expected to be large inside the plunging region 
where the field is forced to become mostly radial or above the disk within 
the corona where the field is mostly toroidal. 



the mass flux, M(r^, t), nominal efficiency, e(rn, /), specific angular 
momentum, j(rii,t), normalized absolute magnetic flux, i>r(rn,t), 
(normalized using the unperturbed initial total flux), and specific 
magnetic flux, T(rH, t), ah measured at the event horizon (r = r^). 
These fluxes have been integrated over the entire range of 9 from 
to n. The quantities M, e and j appear to saturate already at 
t ~ 7000M. However, the magnetic field parameters saturate only 
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Figure 3. Magnetic field lines (red vectors) and magnetic energy density 
(greyscale map) are shown for the fiducial model (A0HR07). Panel (a): 
Snapshot of the magnetic field structure at time 27200M shows that the disk 
is highly turbulent for r > nsco = 6M and laminar for r < nsco- Panel (b): 
Time-averaged magnetic field in the saturated state shows that for r < 9M, 
viz., the region of the flow that we expect to have achieved inflow equilib- 
rium, the geometry of the time-averaged magnetic field closely resembles 
that of a split-monopole. The dashed, vertical line marks the position of the 
ISCO. 

at ~ 12500M. We consider the steady-state period of the disk to 
begin only after all these quantities reach their saturated values. 

The mass accretion rate is quite variable, with root-mean- 
square (rms) fluctuations of order two. The nominal efficiency e 
is fairly close to the NT efficiency, while the specific angular mo- 
mentum J is clearly below the NT value. The results indicate that 
torques are present within the ISCO, but do not dissipate much en- 
ergy or cause significant energy to be transported out of the ISCO. 
The absolute magnetic flux per unit initial absolute flux, 0„ thread- 
ing the black hole grows to about 1 %, which indicates that the mag- 
netic field strength near the black hole is not just set by the amount 
of magnetic flux in the initial torus. This suggests our results are 
insensitive to the total absolute magnetic flux in the initial torus. 
The specific magnetic flux, Y x 0.86 on average. Magnetic stresses 
are relatively weak since T < 1, which implies the mag netic field 
contr ibutes no more than 7% to deviations from NT in j dOammiel 
1 19991) ; see Appendix IaI During the quasi-steady state period, the 
small deviations from NT in / are correlated in time with the magni- 
tude of T. This is consistent with the fact that the specific magnetic 
flux controls these deviations. Also, notice that is roughly con- 
stant in time while T varies in time. This is clearly because M is 
varying in time and also consistent with the fact that T and M are 
anti-correlated in time. 

5.4 Disk Thickness and Fluxes vs. Radius 

Figure [8] shows the time-averaged disk thickness of the fiducial 
model as a function of radius. Both measures of thickness defined 
in 33. H are shown; they track each other As expected, our primary 
thickness measure, \h/r\, is the smaller of the two. This thickness 
measure varies by a small amount across the disk, but it is gen- 
erally consistent with the following fiducial value, viz., the value 
\h/r\ = 0.064 at r = 2r,sco = 12M. 

Figure|9]shows the behavior of various fluxes versus radius for 
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Figure 4. Flow stream lines (red vectors) and rest-mass density (greyscale 
map) are shown for the fiducial model (A0HR07). Panel (a): Snapshot of 
the velocity structure and rest-mass density at time 27200M clearly shows 
MRI-driven turbulence in the interior of the disk. The rest-mass density ap- 
pears more diffusively distributed than the magnetic energy density shown 
in Figure [3ji- Panel (b): Time-averaged streamHnes and rest-mass density 
show that for r < 9M the velocity field is mostly radial with no indication 
of a steady outflow. Time-averaging smooths out the turbulent fluctuations 
in the velocity. The dashed, vertical Hne marks the position of the ISCO. 
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Figure 5. Flow stream lines are shown for the fiducial model (A0HR07). 
Panel (a): Snapshot of the velocity structure at time 27200M clearly shows 
MRI-driven turbulence in the interior of the disk. Panel (b): Time-averaged 
streamlines show that for r < 9M the velocity field is mostly radial. The 
dashed, vertical line marks the position of the ISCO. 

the full 9 integration range (0 to n). We see that the mass accretion 
rate, M, and the specific angular momentum flux, j, are constant 
up to a radius r ~ 9M. This is exactly the distance out to which 
we expect inflow equilibrium to have been established, given the 
inflow velocity and viscous time scale results discussed in 15.21 
The consistency of these two measurements gives us confidence 
that the simulation has truly achieved steady-state conditions inside 
r = 9M. Equally clearly, and as also expected, the simulation is not 
in steady-state at larger radii. 
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Figure 6. The time-averaged, angle-averaged, rest-mass density-weighted 
3-velocities and viscous timescale in the fiducial model (A0HR07) are com- 
pared with the NT model. Angle-averaging is performed over the disk gas 
lying within ±2\h/r\ of the midplane. Top Panel: The orthonormal radial 
3-velocity (solid line), and the analytical GR estimate given in Eq. |B7| of 
Appendix|B](dashed line). Agreement for r > /'isco between the simulation 
and NT model is found when we set a\h/r\- a; 0.00033. At smaller radii, 
the gas dynamics is no longer determined by viscosity and hence the two 
curves deviate. Middle Panel: Shows the orthonormal azimuthal 3-velocity 
(solid line) and the corresponding Keplerian 3-velocity (dashed line). 

Bottom Panel: The inflow equilibrium time scale fie 2r/Vi- (solid line) 

of the disk gas is compared to the analytical GR thin disk estimate (dashed 
line). At r ~ 9M, we see that ~ 2 X lO^M. Therefore, the simulation 
needs to be run for this time period (which we do) before we can reach 
inflow equilibrium at this radius. 

The second panel in Figure |9] shows that the inward angular 
momentum flux, 7in, agrees reasonably well with the NT predic- 
tion. It falls below the NT curve at large radii, i.e., the gas there 
is sub-Keplerian. This is not surprising since we have included the 
contribution of the corona-wind-jet gas which, being at high lat- 
itude, does not rotate at the Keplerian rate. Other quantities, de- 
scribed below, show a similar elfect due to the corona. At the hori- 
zon, j,„ = 3.286, which is 5% lower than the NT value. This devia- 
tion is larger than that found by SOS. Once again, it is because we 
have included the gas in the corona-wind-jet, whereas SOS did not. 

The third panel in Figure |9] shows that the nominal efficiency 
e at the horizon lies below the NT prediction. This implies that the 
full accretion flow (disk-l-corona-l-wind-l-jet) is radiatively less effi- 
cient than the NT model. However, the overall shape of the curve 
as a function of r is similar to the NT curve. The final panel in 
Figure |9] shows the value of T vs. radius. We see that Y » 0.S6 is 
constant out to r ~ 6M. A value of T ~ 1 would have led to 7% 
deviations from NT in j, and only for T ~ 6.0 would deviations 
become 50% (see Appendix |Aj. The fact that Y ~ 0.S6 < 1 in- 
dicates that electromagnetic stresses are weak and cause less than 
7% deviations from NT in j. Note that one does not expect T to be 
constanQ outside the ISCO where the magnetic field is dissipating 

^ We also find that the ideal MHD invariant related to the "isorotation law" 
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Figure 7. Shows for the fiducial model (A0HR07) the time-dependence at 
the horizon of the mass accretion rate, M (top panel); nominal efficiency, e, 
with dashed line showing the NT value (next panel); accreted specific an- 
gular momentum, j, with dashed line showing the NT value (next panel); 
absolute magnetic flux relative to the initial absolute magnetic flux, <I),- 
(next panel); and dimensionless specific magnetic flux, T (bottom panel). 
All quantities have been integrated over all angles. The mass accretion rate 
varies by factors of up to four during the quasi-steady state phase. The nom- 
inal efficiency is close to, but on average slightly lower than, the NT value. 
This means that ffie net energy loss through photons, winds, and jets is be- 
low ffie radiative efficiency of the NT model. The specific angular momen- 
tum is clearly lower than ffie NT value, which implies ffiat some stresses 
are present inside the ISCO. The absolute magnetic flux at the black hole 
horizon grows until it saturates due to local force-balance. The specific mag- 
netic flux T < 1, indicating that electromagnetic stresses inside the ISCO 
are weak and cause less ffian 7% deviations from NT in j. 

due to MHD turbulence and the gas is forced to be nearly Keplerian 
despite a sheared magnetic field. 

As we have hinted above, we expect large differences between 
the properties of the gas that accretes in the disk proper, close to the 
midplane, and that which flows in the corona-wind-jet region. To 
focus just on the disk gas, we show in Figure[TO]the same fluxes as 
inFigure|9] except that we have restricted the (Grange to n/2±2\h/r\. 
The mass accretion rate is no longer perfectly constant for r < 9M. 
This is simply a consequence of the fact that the flow streamlines 
do not perfectly follow the particular constant 2\h/r\ disk boundary 
we have chosen. The non-constancy of M does not significantly 
affect the other quantities plotted in this figure since they are all 
normalized by the local M. 

The specific angular momentum, specific energy, and specific 
magnetic flux are clearly shifted closer to the NT values when we 
restrict the angular integration range. Compared to the NT value, 
viz., 7nt(''h) = 3.464, the fiducial model gives j(ryi) = 3.363 (2.9% 
less than NT) when integrating over ±2\h/r\ around the midplane 
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Figure 8. The time-averaged scale-height, \h/r\, vs. radius in the fiducial 
model (A0HR07) is shown by the solid lines. The above-equator and below- 
equator values of the disk thickness are \h/r\ ~ 0.04-0.06 in the inflow 
equilibrium region r < 9M. We use the specific value of \h/r\ = 0.064 as 
measured at r = 2risco (light dashed lines) as a representative thickness for 
the entire flow. Twice this representative thickness (thick dashed lines) is 
used to fix the 9 range of integration for averaging when we wish to focus 
only on the gas in the disk instead of the gas in the corona-wind-jet. The 
root mean square thickness (/i/r)rms ~ 0.07-0.13 is shown by the dotted 
lines. 



(i.e., only over the disk gas) and gives j(rn) = 3.266 (5.7% less 
than NT) when integrating over all (i.e., including the corona- 
wind-jet). Even though the mass accretion rate through the corona- 
wind-jet is much lower than in the disk, still this gas contributes 
essentially as much to the deviation of the specific angular momen- 
tum as the disk gas does. In the case of the specific magnetic flux, 
integrating over ±2|A/r| around the midplane we find Y a 0.45, 
while when we integrate over all angles T « 0.86. The iGammid 
( Il999h model of an equatorial (thin) magnetized flow within the 
ISCO shows that deviations in the specific angular momentum are 
determined by the value of T. We find that the measured values of 
T are able to roughly predict the measured deviations from NT in 
J- 

In summary, a comparison of Figure |9] and Figure [TOl shows 
that all aspects of the accretion flow in the fiducial simulation agree 
much better with the NT prediction when we restrict our attention 
to regions close to the midplane. In other words, the gas in the disk 
proper, defined here as the region lying within ±2\hlr\ of the mid- 
plane, is well described by the NT model. The deviation of the an- 
gular momentum flux ji„ or j at the horizon relative to NT is < 3%, 
similar to the deviation found by SOfl while the nominal efficiency 
e agrees to within ~ 1%. 



of field fines, (r) = (/ / d0d<p -f^WB't' - f^Bn) / (/ / dBdij, VEMlll 
is Ke plerian outside the ISCO and is (as predicted by the iGammij 
T999I model) roughly constant from ffie ISCO to the ho rizon (see also 
McKinnev & Gammiel2004l ; lMcKinnev & Narav"ai]|2007 j) . 



The quantities j\„ and j are nearly equal at the horizon in the calculations 
reported here whereas they were dift'erent in SOS. This is because SOS used 
an alternate definition of If we had used that definition here, we would 
have found a deviation of ~ 2% in just as in SOS 
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Figure 9. Mass accretion rate and specific fluxes are shown as a function 
of radius for the fiducial model (A0HR07). From top to bottom the pan- 
els show: Top Panel: mass accretion rate; Second Panel: the accreted spe- 
cific angular momentum, j (dotted line), yin (solid line), and the NT profile 
(dashed line); Third Panel: the nominal efficiency e (solid line) and the NT 
profile (dashed line); Bottom Panel: the specific magnetic flux T. For all 
quantities the integration range includes all 6. The mass accretion rate and 
] are roughly constant out to r ~ 9M, as we would expect for inflow equi- 
librium. The profile of ]\a lies below the NT value at lai'ge radii because 
we include gas in the slowly rotating corona. At the horizon, j and e aie 
modestly below the corresponding NT values. The quantity T ~ 0.86 and 
is roughly constant out to r ~ 6M, indicating that electromagnetic stresses 
are weak inside the ISCO. 
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Figure 10. Similar to Figure |9] but here the integration range only includes 
angles within ±2\hlr\ = ±0.128 radians of the midplane. This allows us to 
focus on the disk gas. The mass accretion rate is no longer constant because 
streamlines are not precisely radial. The quantities shown in the second and 
third panels are not aft'ected by the non-constancy of M because they are 
ratios of time-averaged fluxes within the equatorial region and are related 
to ideal MHD invariants. As compared to Figure |9] here we find that j, 
Jin, and e closely follow the NT model. For example, jifw) = 3.363 is 
only 2.9% less than NT. This indicates that the disk and coronal regions 
behave quite differently. As one might expect, the disk region behaves like 
the NT model, while the corona-wind-jet does not. The specific magnetic 
flux is even smaller than in Figure |2]and is T ~ 0.45, which indicates that 
electromagnetic stresses ai'e quite weak inside the disk near the midplane. 



5.5 Comparison with Gammie (1999) Model 

Fig ure 1111 shows a comparison between the fiducial model and 
the iGammid j 19991) model of a magnetized thin accretion flow 
within the ISCO (see also Appendix |AI. Quantities have been in- 
tegrated within ±2\h/r\ of the midplane and time-averaged over a 
short period from t = 17400M to t = 18400M. Note that time- 
averaging b^, po, etc. over long periods can lead to no consistent 
comparable solution if the value of T varies considerably during 
the period used for averaging. Also, note that the presence of ver- 
tical stratification, seen in Figures [9] and \T0\ showing that T de- 
pends upon height, means the vertical-averaging used to obtain T 
can sometimes make it difficult to compare the simulations with 
thelOammie (1999) model which has no vertical stratification. In 
particular, using equation ([24} over this time period, we find that 
T a 0.2 , 0.3, 0.44, 0.7, 0.8 for integrations around the midplane 
of, respectively, ±0.01, ± 0.05, ± 2\h/r\, ± 7t/4, ± 7t/2, with 
best matches to the Gammie model (i.e. and other quantities 
match) using an actual value of T = 0.2, 0.33, 0.47, 0.8, 0.92. 
This indicates that stratification likely causes our diagnostic to un- 
derestimate the best match with the Gammie model once the inte- 
gration is performed over highly-stratified regions. However, the 
consistency is fairly good considering how much T varies with 
height. 

Overall, Figure [TT] shows how electromagnetic stresses con- 
trol the deviations from NT within the ISCO. The panels with D[j] 



and D[e] show how the electromagnetic flux starts out large at the 
ISCO and drops to nearly zero on the horizon. This indicates the 
electromagnetic flux has been converted into particle flux within 
the ISCO by ideal (non-dissipative) electromagnetic stresses^. The 
simulated magnetized thin disk agrees quite well with the Gam- 
mie solutio n, in contrast to the relatively poor agreement found for 
thick disks jMcKinnev & Gammij|2004l) . Only the single parame- 
ter T determines the Gammie solution, so the agreement with the 
value and radial dependence among multiple independent terms is 
a strong validation that the Gammie model is working well. Nev- 
ertheless, there are some residual deviations near the ISCO where 
the thermal pressure dominates the magnetic pressure. Even if de- 
viations from NT are present right at the ISCO, the total deviation 
of the particle flux betw een the ISCO an d horizon equals the de- 
viation predicted by thelGammid ( 1 19991) model, as also found in 
iMcKinnev & Gammid ( |2004|) for thick disks. This indicates that 



Gammid d 19991) ^ model accurately predicts the effects of elec- 



the 

tromagnetic stresses inward of the ISCO. 

Finally, note that the electromagnetic stresses within the ISCO 
are ideal and non-dissipative in the Gammie model. Since the flow 
within the ISCO in the simulation is mostly laminar leading to weak 
non-ideal (resistive or viscous) effects, the dissipative contribution 



' This behavior is just like that seen in ideal MHD jet solutions, but inverted 
with radius. 
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Figure 11. Comparison between the accretion flow (within ±2\h/r\ around 
the midplane) in the fiducial model (A0HR07), shown by solid lines, and 
the mod el of a magnetize d thin accretion disk (inflow solution) within the 
ISCO bv lGamrni3 l ll999l) . shown by dashed lines. In all cases the red ver- 
tical line shows the location of the ISCO. Top-left panel: Shows the radial 
4-velocity, where the Gammie solution assumes = at the ISCO. Fi- 
nite thermal effects lead to non-zero li at the ISCO for the simulated disk. 
Bottom-left panel: Shows the rest-mass density (po, black line), the inter- 
nal energy density (lig, magenta line), and magnetic energy density (h^ 12, 
green line). Top-right and bottom-right panels: Show the percent deviations 
from NT for the simulations and Gammie solution for the specific particle 
kinetic flux (m^, black line), specific enthalpy flux ({Ug + Pg)"^/po, magenta 
line), and specific electromagnetic flux ((b^ u'' Uf^ - b'' bi^,) / (pgu''), green line), 
where for j we use = (p and for e we use fi = t. As usual, the simu- 
lation result for the specific fluxes is obtained by a ratio of flux integrals 
instead of the direct ratio of flux densities. The total specific flux is con- 
stant vs. radius and is a sum of the particle, enthalpy, and electromagnetic 
terms. This figure is comparable to Fig. 10 for a thick (\h/r\ ~ 0.2-0.25) 
disk in iMcKinnev & Gammid fiooj) . Finite thermal pressure effects cause 
the fiducial model to deviate from the inflow solution near the ISCO, but 
the solutions rapidly converge inside the ISCO and the differences between 
the simulation result and the Gammie model (relative to the total specific 
angular momentum or energy) are less than 0.5%. 

(which could lead to radiation) can be quite small. An exception 
to this is the presence of extended current sheets, present near the 
equator within the ISCO in the simulations, whose dissipation re- 
quires a model of the (as of yet, poorly understood) rate of rela- 
tivistic reconnection. 

5.6 Luminosity vs. Radius 

Figure [T2I shows radial profiles of two measures of the disk lumi- 
nosity: L(< r)/M, which is the cumulative luminosity inside radius 
r, and d(L/M)/dln r, which gives the local luminosity at r. We see 
that the profiles from the simulation are quite close to the NT pre- 
diction, especially in the steady-state region. As a way of measur- 
ing the deviation of the simulation results from the NT model, we 
estimate what fraction of the disk luminosity is emitted inside the 
ISCO; recall that the NT model predicts zero luminosity here. The 
fiducial simulation gives L(< nsco)/^ = 0.0021, which is 3.5% 



5 10 15 20 

r/M 

Figure 12. Luminosity per unit rest-mass accretion rate vs. radius (top 
panel) and the logarithmic derivative of this quantity (bottom panel) are 
shown for the fiducial model (A0HR07). The integration includes all 9 an- 
gles. The simulation result (solid lines, truncated into dotted lines outside 
the radius of inflow equilibrium) shows that the accretion flow emits more 
radiation than the NT prediction (dashed lines) at small radii. However, the 
excess luminosity within the ISCO is only Lin ~ 3.5%, where ?[NT] is the 
NT efficiency at the horizon (or equivalently at the ISCO). 

of the nominal efficiency e[NT] = 0.058 of a thin NT disk around 
a non-spinning black hole. This shows that the excess luminosity 
radiated within the ISCO is quite small. The relative luminosity 
within the ISCO is = 3.5% and the relative luminosity within 
the inflow equilibrium region is L(,„( = 8.0%. Hence, we conclude 
that, for accretion disks which are as thin as our fiducial model, 
viz., \h/r\ ~ 0.07, the NT model provides a good description of the 
luminosity profile. 

5.7 Luminosity from Disk vs. Corona- Wind- Jet 

The fiducial model described so far includes a tapering of the cool- 
ing rate as a function of height above the midplane, given by the 
function S [9] (see equation|5]l. We introduced this taper in order to 
only cool bound (-u,(po + Ug + pg + b-)/po < 1) gas and to avoid 
including the emission from the part of the corona- wind-jet that is 
prone to excessive numerical dissipation due to the low resolution 
used high above the accretion disk. This is a common approach that 
others have also taken when performing GRMHD simulations of 
thin disks (N09, NIO). However, since our tapering function does 
not explicitly refer to how bound the gas is, we need to check that 
it is consistent with cooling only bound gas. We have explored this 
question by re-running the fiducial model with all parameters the 
same except that we turned off the tapering function altogether, i.e., 
we set S[6] = 1 . This is the only model for which the tapering func- 
tion is turned off. 

Figure [13] shows a number of luminosity profiles for the fidu- 
cial model and the no-tapering model. This comparison shows that, 
whether or not we include a taper, the results for the luminosity 
from all the bound gas is nearly the same. Without a tapering, there 
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is some luminosity at high latitudes above ±8|/i/r| corresponding 
to emission from the low-density jet region (black solid line). This 
region is unbound and numerically inaccurate, and it is properly 
excluded when we use the tapering function. Another conclusion 
from the above test is that, as far as the luminosity is concerned, 
it does not matter much whether we focus on the midplane gas 
((n/2) ± 2\h/r\) or include all the bound gas. The deviations of the 
luminosity from NT in the two cases are similar - changes in the 
deviation are less than 1%. 

An important question to ask is whether the excess luminosity 
from within the ISCO is correlated with, e.g., deviations from NT 
in J, since D[j] could then be used as a proxy for the excess lu- 
minosity. We investigate this in the context of the simulation with 
no tapering. For an integration over ±2\h/r\ around the midplane 
(which we identify with the disk component), or over all bound 
gas, or over all the gas (bound and unbound), the excess luminos- 
ity inside the ISCO is !,„ = 3.3%, 4.4%, 5.4%, and the deviation 
from NT in ; is D[j] = -3.6%, - 6.7%, - 6.7%, respectively. We 
ignore the luminosity from unbound gas since this is mostly due 
to material in a very low density region of the simulation where 
thermodynamics is not evolved accurately. Considering the rest of 
the results, we see that D[j] is 100% larger when we include bound 
gas outside the disk compared to when we consider only the disk 
gas, whereas the excess luminosity increases by only 32%. There- 
fore, when we compute j by integrating over all bound gas and 
then assess the deviation of the simulated accretion flow from the 
NT model, we strongly overestimate the excess luminosity of the 
bound gas relative to NT. A better proxy for the latter is the devi- 
ations from NT in j integrated only over the disk component (i.e. 
over ±2\h/r\ around the midplane). 

Furthermore, we note that the gas that lies beyond ±2\h/r\ 
from the disk midplane consists of coronal gas, which is expected 
to be optically thin and to emit a power-law spectrum of photons. 
For many applications, we are not interested in this component but 
rather care only about the thermal blackbody-like emission from 
the optically-thick region of the disk. For such studies, the most ap- 
propriate diagnostic from the simulations is the radiation emitted 
within ±2\h/r\ of the midplane. According to this diagnostic, the 
excess emission inside the ISCO is only Z,j„ = 3.4% in the model 
without tapering, and 3.5% in the fiducial model that includes ta- 
pering. 

Lastly, we consider variations in the cooling timescale, TcooI , 
which is another free parameter of our cooling model that we gen- 
erally set to 2n/Qi<_- However, we consider one model that is oth- 
erwise identical to the fiducial model except we set Tcooi to be five 
times shorter so that the cooling rate is five times faster. We find 
that Li„ = 4.2%, which is slightly larger than the fiducial model 
with Li„ = 3.5%. Even though the cooling rate is five times faster 
than an orbital rate, there is only 20% more luminosity from within 
the ISCO. This is likely due to the flow within the ISCO being 
mostly laminar with little remaining turbulence to drive dissipation 
and radiation. 



6 CONVERGENCE WITH RESOLUTION AND BOX SIZE 

The fiducial model described earlier was computed with a numeri- 
cal resolution of 256 x 64 x 32, using an azimuthal wedge of ;r/2. 
This is to be compared with the simulation described in SOS, which 
made use of a 512 x 128 x 32 grid and used a 7r/4 wedge. These 
two runs give very similar results, suggesting that the details of the 
resolution and wedge size are not very important. To confirm this. 



T r 



8x10-3 - 




Figure 13. Shows enclosed luminosity vs. radius for models with different 
cooling prescriptions and integration ranges. The black dashed line cor- 
responds to the NT model. The luminosity for the fiducial model A0HR07, 
which includes a tapering of the cooling with disk height as described in ^ 
is shown integrated over ±2\h/r\ from the midplane (black dotted hne), in- 
tegrated over all bound gas (black long dashed line), and integrated over all 
fluid (black solid line). Essentially all the gas is bound and so the black solid 
and long dashed lines are indistinguishable. The red lines are for a model 
that is identical to the fiducial mn, except that no tapering is applied to the 
cooling. For this model the lines ai'e: red solid line: all angles, all fluid; red 
dotted line: ±2\h/r\ ai'ound the midplane; red long dashed line: all bound 
gas. The main result is that the luminosity from bound gas is nearly the 
same (especially at the ISCO) whether or not we include tapering (compare 
the red long dashed line and the black long dashed line). 



we have run a number of simulations with different resolutions and 
wedge angles. The complete list of runs is: 256 x 64 x 32 with 
A(f> = ;r/2 (fiducial run, model A0HR07)), 512 x 128 x 32 with 
A(f> = 7r/4 (SOS, model CO), 256 x 12S x 32 with A(p = n/2 (model 
C6), 256 X 32 X 32 with A(f> = 7t/2 (model C5), 256 x 64 x 64 with 
A(f) = n/2 (model C4), 256 x 64 x 64 with A(/> = n (model C2), 
256 X 64 X 16 with = ;r/2 (model C3), and 256 x 64 x 16 with 
A(f> = n/4 (model CI). 

Figure [T4l shows the accreted specific angular momentum, j, 
ingoing component of the specific angular momentum, ji„, and the 
nominal efliciency e as functions of radius for all the models used 
for convergence testing. Figure [15] similarly shows the cumulative 
luminosity L(< r)/M and differential luminosity d(L/ M)/dlnr as 
functions of radius. The overwhelming impression from these plots 
is that the sequence of convergence simulations agree with one an- 
other quite well. Also, the average of all the runs matches the NT 
model very well; this is especially true for the steady-state region of 
the flow, r < 9M. Thus, qualitatively, we conclude that our results 
are well-converged. 

For more quantitative comparison. Figure [16] shows the pro- 
file of / vs r for the various models, this time with each model 
separately identified. It is clear that / has converged, since there are 
very minor deviations from our highest resolution/largest box size 
to our next highest resolution/next largest box size. All other quan- 
tities, including e, ji„, and T are similarly converged. The model 
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Table 2. Convergence Study 



Model Name \h/r\ M e D[e] j D[j] D[j,^] s Li„ lOOd^ T 



±2\h/r\ 




























A0HR07 


0.064 


0.066 


0.058 


-0.829 


3.363 


-2.913 


3.355 


-3.153 


3.363 


0.035 


0.080 


1.355 


0.450 


CO 


0.052 


0.064 


0.059 


-2.708 


3.363 


-2.905 


3.351 


-3.259 


3.363 


0.019 


0.042 


0.399 


0.359 


CI 


0.063 


0.041 


0.056 


2.204 


3.415 


-1.406 


3.408 


-1.621 


3.415 


0.034 


0.072 


0.584 


0.223 


C2 


0.062 


0.063 


0.058 


-0.699 


3.360 


-3.016 


3.333 


-3.789 


3.360 


0.047 


0.109 


0.925 


0.323 


C3 


0.061 


0.026 


0.058 


-1.969 


3.358 


-3.060 


3.339 


-3.610 


3.358 


0.039 


0.087 


0.727 


0.449 


C4 


0.061 


0.054 


0.058 


-1.406 


3.385 


-2.286 


3.378 


-2.489 


3.385 


0.019 


0.018 


0.714 


0.296 


C5 


0.052 


0.008 


0.055 


3.955 


3.417 


-1.364 


3.427 


-1.067 


3.417 


0.034 


0.070 


0.315 


0.322 


C6 


0.065 


0.088 


0.059 


-3.355 


3.355 


-3.155 


3.333 


-3.778 


3.355 


0.054 


0.103 


0.933 


0.256 


Aiie 




























A0HR07 


0.064 


0.074 


0.054 


4.723 


3.266 


-5.717 


3.281 


-5.275 


3.266 


0.035 


0.053 


6.677 


0.863 


CO 


0.052 


0.071 


0.057 


0.738 


3.312 


-4.392 


3.291 


-5.002 


3.312 


0.032 


0.049 


2.18 


0.480 


CI 


0.063 


0.042 


0.055 


3.680 


3.398 


-1.894 


3.389 


-2.178 


3.398 


0.037 


0.062 


0.940 


0.142 


C2 


0.062 


0.069 


0.056 


1.664 


3.307 


-4.541 


3.262 


-5.846 


3.307 


0.053 


0.080 


5.757 


0.710 


C3 


0.061 


0.029 


0.057 


0.893 


3.325 


-4.008 


3.305 


-4.580 


3.325 


0.042 


0.064 


1.401 


0.309 


C4 


0.061 


0.057 


0.057 


0.811 


3.358 


-3.075 


3.351 


-3.255 


3.358 


0.020 


0.009 


2.690 


0.359 


C5 


0.052 


0.009 


0.053 


7.687 


3.338 


-3.636 


3.353 


-3.218 


3.338 


0.039 


0.059 


0.726 


0.243 


C6 


0.065 


0.092 


0.058 


-1.813 


3.334 


-3.761 


3.306 


-4.560 


3.334 


0.057 


0.086 


5.091 


0.534 




Figure 14. This plot shows j, and e for a sequence of simulations that 
are similar to the fiducial run (A0HR07), viz., \h/r\ « 0.07, a/M = 0, but 
use different radial resolutions, or resolutions, or box sizes. The integra- 
tion range in 8 is over ±2\h/r\ around the midplane. Only the region of the 
flow in inflow equiUbrium, 2M < r < 9M, is shown in the case of j. The 
different lines are as follows: black dashed line: NT model; black solid hne: 
fiducial model A0HR07; blue solid line: model CO (S08); magenta dotted 
line: model CI; magenta solid line: model C2; red dotted line: model C3; 
red solid line: model C4; green dotted line: model C5; green solid line: 
model C6. Note that changes in the numerical resolution or other computa- 
tional parameters lead to negligible changes in the values of j, and e in 
the region of the flow that is in inflow equilibrium, r < 9M. For r 5; 9M, 
the flow has not achieved steady state, which explains the large deviations 
in e. Only the lowest resolution models are outliers. 



Figure 15. Similar to Figure [14] but for the normalized luminosity, L(< 
r)/M, and its logarithmic derivative, d(L/M)/dln r, both shown vs. radius. 
We see that all the models used to test convergence show consistent lumi- 
nosity profiles over the region that is in inflow equilibrium, r < 9M. The 
well-converged models have Li„ <, 4%, which indicates only a low level of 
luminosity inside the ISCO. 

with A'^ = 64 shows slightly less deviations from NT in j than our 
other models. However, it also shows slightly higher luminosity 
than our other models. This behavior is likely due to the stochastic 
temporal behavior of all quantities vs. time, but this could also be 
due to the higher ^-resolution causing a weaker ordered magnetic 
field to be present leading to weaker ideal electromagnetic stresses, 
smaller deviations from NT in j within the ISCO, but with the re- 
maining turbulent field being dissipated giving a higher luminosity. 
The A'^ resolution appears to be the limit on our accuracy. 
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Figure 16. This is a more detailed version of Figure [T4l showing y vs r for 
individually labeled models. The models correspond to the fiducial reso- 
lution (solid lines), a higher resolution ntn (dot-dashed Unes), and a lower 
resolution run (dotted lines). Generally, there are only minor differences 
between the fiducial and higher resolution models. 



Further quantitative details are given in Table|2l where we list 
numerical results for all the convergence test models, with the 6 
integration performed over both ±2|/i/r| around the midplane and 
over all angles. We see that there are some trends as a function of 
resolution and/or A0. Having only 32 cells in 9 or 16 cells in gives 
somewhat poor results, so these runs are under-resolved. How- 
ever, even for these runs, the differences are not large. Note that 
T reaches a steady-state much later than all other quantities, and 
our C? (where ? is through 6) models did not run as long as the 
fiducial model. This explains why Y is a bit lower for the C? mod- 
els. Overall, we conclude that our choice of resolution 256x 64x 32 
for the fiducial run (A0HR07) is adequate to reach convergence. 



7 DEPENDENCE ON BLACK HOLE SPIN AND DISK 
THICKNESS 

In addition to the fiducial model and the convergence runs de- 
scribed in previous sections, we have run a number of other sim- 
ulations to explore the effect of the black hole spin parameter a/M 
and the disk thickness \hlr\ on our various diagnostics: j, jj„, e, the 
luminosity, and Y. We consider four values of the black hole spin 
parameter, viz., a/M = 0, 0.7, 0.9, 0.98, and four disk thicknesses, 
viz., = 0.07, 0.1, 0.2, 0.3. We summarize here our results for 
this 4x4 grid of models. 

Geometrically thick disks are expected on quite general 
grounds to deviate from the standard thin disk model. The inner 
edge of the disk, as measured for instance by the location of the 
sonic point, is expected to deviate from the ISCO, the shift scal- 
ing roughly as |ri„ - nscol °^ (Cs/vk)^ (c\ is sound speed, where 
cl = Ypglipa + Ug + p g)). This effect i s seen in h ydrodynamic mod- 



els of thick disks, e.g 
( I2OIOI) , where it is shown that n, 



Narayan et al.l l ll997h and lAbramowicz et al.l 



can move either inside or outside 



the ISCO; it moves inside when a is small and outside when a is 
large. In either case, these hydrodynamic models clearly show that, 
as \hlr\ — > 0, i .e., as c^/dk — » 0, th e solution always tends toward 
the NT model jShafee et alj2008bh . 

While the hydrodynamic studies mentioned above have driven 
much of our intuition on the behavior of thick and thin disks, it is an 
open question whether or not the magnetic field plays a significant 
role. In principle, magnetic effects may cause the solution to deviate 
significantly from the NT model even in the limit \hlr\ —> ( KrolikI 
ll999l : lGammiJl999l) . One of the major goals of the present paper is 
to investigate this question. We show in this section that, as \h/r\ — > 
0, magnetized disks do tend toward the NT model. This statement 
appears to be true for a range of black hole spins. We also show 
that the specific magnetic flux Y inside the ISCO decreases with 
decreasing \h/r\ and remains quite small. This explains why the 
magnetic field does not cause significant deviations from NT in 
thin disks. 

Figure [T7] shows the specific angular momentum, j, and the 
ingoing component of this quantity, vs. radius for the 4x4 grid 
of models. The 6 integral has been taken over ±2\h/r\ around the 
midplane in order to focus on the equatorial disk properties. The 
value of J is roughly constant out to a radius well outside the ISCO, 
indicating that we have converged solutions in inflow equilibrium 
extending over a useful range of radius. As discussed in section [33] 
inflow equilibrium is expected within r/M = 9, 7, 5.5, 5, respec- 
tively, for a/M = 0, 0.7, 0.9, 0.98. This is roughly consistent with 
the radius out to which the quantity j (integrated over all angles) 
is constant, and this motivates why in all such plots we only show 
J over the region where the flow is in inflow equilibrium. The four 
panels in Figure[T7]show a clear trend, viz., deviations from NT are 
larger for thicker disks, as expected. Interestingly, for higher black 
hole spins, the relative deviations from NT actually decrease. 

Figure [TS] shows the nominal efficiency, e, as a function of 
radius for the 4x4 grid of models. Our thickest disk models (\h/r\ x 
0.3) do not include cooling, so the efficiency shown is only due to 
losses by a wind-jet. We see that the efficiency is fairly close to 
the NT value for all four thin disk simulations with \h/r\ ~ 0.07; 
even in the worst case, viz., a/M = 0.98, the deviation from NT is 
only ~ 5%. In the case of thicker disks, the efficiency shows larger 
deviations from NT and the profile as a function of radius also looks 
different. For models with \h/r\ x 0.3, there is no cooling so large 
deviations are expected. 

Figure [T9I shows the luminosity, L(< r)/ M, vs. radius for our 
4x4 grid of models, focusing just on the region that has reached 
inflow equilibrium. The luminosity is estimated by integrating over 
all 6 angles. Our thickest disk models {\h/r\ « 0.3) do not include 
cooling and so are not plotted. The various panels show that, as 
\h/r\ — » 0, the luminosity becomes progressively closer to the NT 
result in the steady state region of the flow near and inside the 
ISCO. Thus, once again, we conclude that the NT luminosity pro- 
file is valid for geometrically thin disks even when the accreting 
gas is magnetized. 

A figure (not shown) that is similar to Figure [T7] but for the 
specific magnetic flux indicates that Y < 1 within ±2\h/r\ near 
the ISCO for all black hole spins and disk thic knesses. For our 
thinnest models, Y < 0.45, for which the model o f|G amrm 3l ll999t) 
predicts that the specific angular momentum will deviate from 
NT by less than 1.9%, 3.0%, 3.8%, 4.2% for black hole spins 
a/M = 0, 0.7, 0.9, 0.98, respectively (see Appendix |A). The nu- 
merical results from the simulations show deviations from NT that 
are similar to these values. Thus, overall, our results indicate that 
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Figure 17. The net accreted specific angular momentum, j (the nearly 
liorizontal dotted lines), and tlie ingoing component of this quantity, yi„ 
(the sloping curved lines), as a function of radius for the 4x4 grid 
of models. Each panel corresponds to a single black hole spin, a/M = 
0, 0.7, 0.9, or 0.98, and shows models with four disk thicknesses, \h/r\ « 
0.07, 0.1, 0.2, 0.3 (see legend). The 6 integral has been taken over ±2\h/r\ 
around the midplane. In each panel, the thin dashed black line, marked by 
two circles which indicate the location of the horizon and the ISCO, shows 
the NT solution for As expected, we see that thicker disks exhibit larger 
deviations froin NT. However, as a function of spin, there is no indication 
that deviations froin NT become any larger for larger spins. In the case of 
the thinnest models with \h/r\ x 0.07, the NT model works well for gas 
close to the midplane for all spins. 



electromagnetic stresses are weak inside the ISCO for geometri- 
cally thin disks. 

Finally, for all models we look at plots (not shown) of M{< r) 
(mass enclosed within radius), M{r) (total mass accretion rate vs. 
radius), and [/j/r](r) (disk scale-height vs. radius). We find that 
these are consistently flat to the same degree and to the same ra- 
dius as the quantity ;(r) is constant as shown in Figure [17] This 
further indicates that our models are in inflow equilibrium out to 
the expected radius. 



7.1 Scaling Laws vs. a/M and 

We now consider how the magnitude of j, e, L(< nsco)> and T 
scale with disk thickness and black hole spin. Table[3]lists numer- 
ical results corresponding to Q integrations over ±2\h/r\ around the 
midplane and over all angle j3 Figure [20] shows selected results 
corresponding to models with a non-rotating black hole for quanti- 
ties integrated over ±2\h/r\. We see that the deviations of various di- 
agnostics from the NT values scale roughly as \h/r\. In general, the 
deviations are quite small for the thinnest model with \h/r\ x 0.07. 



Some thicker disk models without cooling show small or slightly neg- 
ative efficiencies, e, which signifies the accretion of weakly unbound gas. 
This can occur when a magnetic field is inserted into a weakly bound gas in 
hydrostatic equilibrium. 



Figure 18. Similar to Figure [TtI but for the nominal efficiency, e. For thin 
{\hlr\ < 0.1) disks, the results are close to NT for all black hole spins. 
As expected, the thicker models deviate significantly from NT. In part this 
is because the ad hoc cooling function we use in the simulations is less 
accurate for thick disks, and in part because the models with \hlr\ x 0.3 
have no coohng and start with marginally bound/unbound gas that implies 
e ~ 0. The a/M = 0.98 models show erratic behavior at large radii where 
the flow has not achieved inflow equilibrium. 
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Figure 19. Similar to Figure [TT] but for the normalized luminosity, L(< 
r)/M. For thin disks, the luminosity deviates only slightly from NT near 
and inside the ISCO. There is no strong evidence for any dependence on the 
black hole spin. The region at large radii has not reached inflow equilibiium 
and is not shown. 
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Table 3. Black Hole Spin and Disk Thickness Study 



Model Name \hlr\ M e Z)[?] ; D[j} D[;i„] s !,„ Leq lOOj),- T 



±2|/i/r| 




























A0HR07 


0.064 


0.066 


0.058 


-0.829 


3.363 


-2.913 


3.355 


-3.153 


3.363 


0.035 


0.080 


1.355 


0.450 


A7HR07 


0.065 


0.050 


0.107 


-2.899 


2.471 


-4.465 


2.527 


-2.294 


1.220 


0.048 


0.084 


0.919 


0.393 


A9HR07 


0.054 


0.045 


0.156 


-0.157 


2.042 


-2.762 


2.074 


-1.213 


0.523 


0.041 


0.082 


0.455 


0.218 


A98HR07 


0.059 


0.013 


0.225 


3.897 


1.643 


-2.335 


1.679 


-0.199 


0.124 


0.069 


0.127 


0.276 


0.228 


AOHRl 


0.12 


4.973 


0.056 


1.470 


3.138 


-9.424 


3.162 


-8.724 


3.138 


0.084 


0.134 


2.976 


0.871 


A7HR1 


0.09 


2.443 


0.099 


4.808 


2.446 


-5.447 


2.524 


-2.409 


1.184 


0.060 


0.108 


0.909 


0.406 


A9HR1 


0.13 


2.133 


0.142 


9.014 


1.947 


-7.261 


2.124 


1.157 


0.402 


0.068 


0.107 


1.064 


0.466 


A98HR1 


0.099 


2.372 


0.213 


8.810 


1.626 


-3.393 


1.703 


1.200 


0.084 


0.062 


0.112 


0.451 


0.254 


A0HR2 


0.18 


48.286 


0.055 


4.134 


2.774 


-19.916 


2.872 


-17.098 


2.774 


0.167 


0.235 


2.518 


1.235 


A7HR2 


0.16 


31.665 


0.049 


52.330 


2.412 


-6.736 


2.576 


-0.425 


1.081 


0.050 


0.034 


0.919 


0.63 1 


A9HR2 


0.21 


40.603 


0.090 


41.922 


1.946 


-7.315 


2.155 


2.624 


0.309 


0.011 


-0.026 


0.795 


0.557 


A98HR2 


0.18 


29.410 


0.191 


18.496 


1.588 


-5.650 


1.870 


11.117 


0.001 


0.052 


0.068 


0.651 


0.459 


A0HR3 


0.350 


44.066 


-0.003 


104.582 


2.309 


-33.331 


2.408 


-30.473 


2.309 


0.000 


-0.049 


3.039 


1.182 


A7HR3 


0.34 


41.045 


-0.007 


106.384 


1.967 


-23.970 


2.236 


-13.549 


0.557 


0.000 


-0.060 


1.956 


0.746 


A9HR3 


0.341 


35.852 


-0.001 


100.582 


1.722 


-17.999 


1.998 


-4.832 


-0.080 


0.000 


-0.073 


1.437 


0.543 


A98HR3 
All 6 
A0HR07 


0.307 


24.486 


-0.015 


106.206 


1.783 


5.987 


1.886 


12.102 


-0.205 


0.000 


-0.104 


0.369 


0.246 


0.064 


0.074 


0.054 


4.723 


3.266 


-5.717 


3.281 


-5.275 


3.266 


0.035 


0.053 


6.677 


0.863 


A7HR07 


0.065 


0.065 


0.093 


10.132 


2.282 


-11.776 


2.511 


-2.933 


1.012 


0.040 


0.040 


8.496 


1.156 


A9HR07 


0.054 


0.060 


0.135 


13.294 


1.860 


-11.404 


2.209 


5.222 


0.303 


0.035 


0.031 


8.945 


1.299 


A98HR07 


0.059 


0.021 


0.171 


26.735 


1.559 


-7.348 


1.799 


6.905 


-0.065 


0.048 


0.039 


2.460 


0.626 


AOHRl 


0.12 


6.036 


0.053 


6.579 


2.908 


-16.046 


2.980 


-13.961 


2.908 


0.087 


0.110 


8.880 


1.247 


A7HR1 


0.09 


2.907 


0.093 


10.343 


2.344 


-9.376 


2.457 


-4.988 


1.074 


0.068 


0.093 


2.295 


0.525 


A9HR1 


0.13 


2.777 


0.128 


17.577 


1.823 


-13.164 


2.069 


-1.449 


0.254 


0.066 


0.075 


4.256 


0.735 


A98HR1 


0.099 


3.235 


0.197 


15.950 


1.425 


-15.291 


1.880 


11.744 


-0.149 


0.078 


0.094 


6.599 


1.349 


A0HR2 


0.18 


59.025 


0.050 


12.631 


2.465 


-28.830 


2.596 


-25.067 


2.465 


0.164 


0.197 


5.798 


1.771 


A7HR2 


0.16 


41.327 


0.046 


55.244 


2.186 


-15.499 


2.394 


-7.453 


0.851 


0.045 


0.015 


2.679 


0.923 


A9HR2 


0.21 


53.746 


0.085 


45.564 


1.739 


-17.164 


2.058 


-2.001 


0.092 


0.012 


-0.031 


3.799 


1.093 


A98HR2 


0.18 


43.815 


0.154 


34.174 


1.411 


-16.120 


1.876 


11.497 


-0.247 


0.045 


0.033 


2.072 


0.887 


A0HR3 


0.350 


49.207 


-0.004 


106.180 


2.128 


-38.572 


2.220 


-35.919 


2.128 


0.000 


-0.049 


4.724 


1.331 


A7HR3 


0.34 


47.146 


-0.007 


107.191 


1.788 


-30.878 


2.065 


-20.166 


0.377 


0.000 


-0.060 


3.433 


0.952 


A9HR3 


0.341 


42.733 


-0.004 


102.370 


1.530 


-27.117 


1.869 


-10.977 


-0.276 


0.000 


-0.073 


2.652 


0.914 


A98HR3 


0.307 


30.293 


-0.014 


105.873 


1.597 


-5.105 


1.655 


-1.624 


-0.390 


0.000 


-0.104 


0.668 


0.273 



Next, we consider fits of our simulation data as a function of 
black hole spin and disk thickness to reveal if, at all, these two pa- 
rameters control how much the flow deviates from NT. In some 
cases we directly fit the simulation results instead of their devia- 
tions from NT, since for thick disks the actual measurement values 
can saturate independent of thickness leading to large non-linear 
deviations from NT. Before making the fits, we ask how quanti- 
ties might scale with a/M and \hlr\. With no disk present, the ro- 
tational sy mmetry forces any scaling to be an even power of black 
hole spin jMcKinnev & GammidllOOi) . However, the presence of 
a rotating disk breaks this symmetry, and any accretion flow prop- 
erties, such as deviations from NT's model, could depend linearly 
upon alM (at least for small spins). This motivates performing a 
linear fit in ajM. Similarly, the thickness relates to a dimensionless 
speed: Cs/dk ~ \hlr\, while there are several different speeds in the 
accretion problem that could force quantities to have an arbitrary 
dependence on |/?/r|. Although, in principle, deviations might scale 
as some power of \hlr\, we assume here a linear scaling oc |/?/r|. 
This choice is driven partly by simplicity and partly by Figure [20l 
which shows that the simulation results agree well with this scaling. 

These rough arguments motivate obtaining explicit scaling 
laws for a quantity's deviations from NT as a function of a/M and 
\hlr\. For all quantities we use the full 4 x 4 set of models, except 
for the luminosity and efficiency we exclude the two thickest mod- 



els in order to focus on the luminosity for thin disks with cooling. 
We perform a linear least squares fit in both ajM and \hlr\, and 
we report the absolute percent difference between the upper 95% 
confidence limit (C+) and the best-fit parameter value (/) given by 
E = 100IC+ - /I/I/I. Note that if £ > 100%, then the best-fit value 
is no different from zero to 95% confidence (such parameter values 
are not reported). After the linear fit is provided, the value of E is 
given for each parameter in order of appearance in the fit. Only the 
statistically significant digits are shown. 

First, we consider how ele ctromagnetic stresses depend upon 
ajM and |/7/r|.|G ammi has shown that the effects of elec- 

tromagnetic stresses are tied to the specific magnetic flux, T, and 
that for T < 1 there are weak electromagnetic stresses causing only 
minor deviations (less than 12% for j across all black hole spins) 
from NT. Let us consider how T should scale with \hlr\, where 
T = ^IlirlMfB'! ( ^j-{rlMYpoif^ in the equatorial plane and is 
assumed to be constant from the ISCO to the horizon. For sim- 
plicity, let us study the case of a rapidly rotating black hole. First, 
consider the boundary conditions near the ISCO provided by the 
disk, where c's/uk ~ \hlr\ and the Keplerian rotation speed reaches 
~ 0.5. This implies ~ Q.5\hlr\. Second, consider the flow 
that connects the ISCO and the horizon. The gas in the disk be- 
yond the ISCO has /? ~ (Cj/ija)^ ~ 10, but reaches /3 ~ 1 in- 
side the ISCO in any GRMHD simulations of turbulent magne- 
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Figure 20. The relative difference between j in tlie simulation and in the NT 
model (top panel), the relative difference between the nominal efficiency, 
e, and the NT value (middle panel), and the luminosity inside the ISCO 
normalized by the net radiative efficiency of the NT model, Ljn (bottom 
panel), where e[NT] has been evaluated at the horizon (equivalently at the 
ISCO). There is a rough linear dependence on \h/r\ for all quantities, where 
a linear fit is shown as a dotted line in each panel. Note that the thicker disk 
models are not expected to behave like NT, and actually have j roughly 
similar across all spins. For \h/r\ x 0.07, the excess luminosity from within 
the ISCO is less than 4% of the total NT efficiency. 



tized disks, which gives that Cj ~ v/^. Thus, va ~ 0.5|/i/r|. Fi- 
nally, notice that T ~ 1.46''/ ^Jpo at the horizon where u'' ~ -1 
and r = M. The Keplerian rotation at the ISCO leads to a mag- 
netic field with orthonormal radial (|6,-| ~ |S' |) and toroidal (|fi^|) 
components with similar values near the ISCO and horizon, giving 
Ifil ~ \B,\ ~ \B^\ ~ \b\ and so Y ~ lA\b\/ Further, the Alfven 
3-speed is = \b\l ^b^ + po + Ug + pg ~ \b\ / -y^oo in any massive 
disk, so that T ~ IAva ~ 0.7|/i/r| for a rapidly rotating black hole. 
Extending these rough arguments to all black hole spins at fixed 
disk thickness also gives that T oc -0.8(fl/M) for fl/M < 0.7. These 
arguments demonstrate three points: 1) T » 1 gives b- Ipo » 1, 
implying a force-free magnetosphere instead of a massive accre- 
tion disk ; 2) Y oc |/i/r| ; and 3) Y oc -(a/M). Since the local 
con dition for the magnetic field e jecting mass is b-/po » 1 (see, 
e.g., lKomissarov & Barkovll2009l) , this shows that Y ~ 1 defines a 
boundary that the disk component of the flow cannot significantly 
pass beyond without eventually incurring disruption by the strong 
magnetic field within the disk. 

We now obtain the actual fit, which for an integration over 
±2\h/r\ gives 



Y a 0.7 + 



, a 
-0.6—, 
M 



(37) 



with E = 33%, 70%, 40%, indicating a reasonable fit. There is 
essentially 100% confidence in the sign of the 1st and 3rd parame- 
ters and 98% confidence in the sign of the 2nd parameter. This fit is 
consistent with our basic analytical estimate for the scaling. Since 
most likely Y < 0.9 in the limit that \h/r\ — > across all black 



hole spins, the electromagnetic stresses are weak and cause less 
than 12% deviation from NT in j. This means that NT solution is 
essentially recovered for magnetized thin disks. For an integration 
over all angles, T x \ with E = 35%, and there is no statistically 
significant trend with disk thickness or black hole spin. The value 
of Y ~ 1 is consistent with the presence of the highly-magnetized 
corona-wind-jet above the disk component jMcKinnev & Gammid 
l2004h . 

Next, we consider whether our simulations can determine the 
equilibrium value of the black hole spin as a function of \h/r\. The 
spin evolves as the black hole accretes mass, energy, and angular 
momentum, and it can stop evolving when these come into a cer- 
tain balance leading to d(a/M)/dt = (see equation [Tsl. In spin- 
equilibrium, the spin-up parameter s = j - 2(a/M)e has .v = 
and solving for a gives the equilibrium spin a^^/M = j/(2e). For 
the NT solution, s is fairly linear for a/M > and a^q/M = 1. In 
appendix [a] we note that for Y ~ 0.2-1 that the deviations from 
NT roughly scale as Y. Since Y oc \h/r\, one expects s to also 
roughly scale with \h/r\. This implies that deviations from NT in 
the spin equilibrium should scale as \h/r\. Hence, one should have 
1 - Ceq/M oc \h/r\. 

Now we obtain the actual fit. We consider two types of fits. In 
one case, we fit s (with fluxes integrated over all angles) and solve 
s = for a^^lM. This gives 



s a 3.2-2.5 



-2.9ii, 



(38) 



with E = 8%, 36%, 8%, indicating quite a good fit. There is an es- 
sentially 100% confidence for the sign of all parameters, indicating 
the presence of well-defined trends. Solving the equation s = Q for 
a/M shows that the spin equilibrium value, a^^/M, is given by 

h 



1 i a -0.08 + O.i 

M 



(39) 



In the other case, we fit j/(2s) and re-solve for a^q/M, which gives 
directly 

h 

a; -0 10 + 04 

M 



(40) 



with E = 9%, 38% with a 99.99% confidence in the sign of the \h/r\ 
term. Both of these procedures give a similar fit (the first fit is sta- 
tistically better) and agree within statistical errors, which indicates 
a linear fit is reasonable. For either fit, one should set a^q/M = 1 
when the above formula gives a^q/M > 1 to be consistent with our 
statistical errors and the correct physics. Note that the overshoot 
Ocq/M > 1 in the fit is consistent with a linear extrapolation of the 
NT dependence of i- for a/M > 0, which also overshoots in the 
same way due to the progressively non-linear behavior of s above 
a/M « 0.95. 

These spin equilibrium fits imply that, within our statistical 
errors, the spin can reach a^q/M — > 1 as \h/r\ — > 0. Thus, our re- 
sults are consistent with NT by allowing maximal black hole spin 
for thin disk^ Our results are also roughly consistent with the 
thick disk 1-loop field geometry study by Gammie et al. (2004j). 
Using our definition of disk thickness, their model had \h/r\ ~ 0.2- 
0.25 and they found a^q/M ~ 0.9, which is roughly consistent with 
our scaling law. The fit is also consistent with results for even 
thicker disks (\h/r ~ . 4 near the h orizon) with a^q/M ~ 0.8 
dAbramowicz et al.|[l978l : |PoDham"&Gammi&,1998() . 



" Here, we do not include black h ole spin chang es by photon capture, 
which gives a limit ofa^q/M = 0.998 jThomell97# . 
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Overall, the precise scaling relations given for T and a^^ 
should be considered as suggestive and preliminary. More work is 
required to test the convergence and generality of the actual coef- 
ficients. While we explicitly tested convergence for the ajM = 
fiducial model, the other ajM were not tested as rigorously. A 
potential issue is that we find the saturated state has fewer cells 
per (vertical magnetic field) fastest growing mode for the axisym- 
metric MRI in models with a/M = 0.9,0.98 than in models with 
fl/M = 0, 0.7 due to the relative weakness of the vertical field in the 
saturated state for the high spin models. However, both the rough 
analytical arguments and the numerical solutions imply that elec- 
tromagnetic stresses scale somewhat linearly with black hole spin. 
This consistency suggests that many measurements for the simula- 
tions, such as T and Oeq, may be independent of smallness of the 
vertical field. This fact could be due to these quantities being only 
directly related to the radial and toroidal magnetic field strengths 
rather than the vertical magnetic field strength. Further, our thick 
disk models resolve the axisymmetric MRI better than the thinnest 
disk model. This suggests that the scaling of T and cieq with disk 
thickness may be a robust result. 

Lastly, we consider how the specific angular momentum, nom- 
inal efficiency, and luminosity from within the ISCO deviate from 
NT as functions of spin and thickness. Overall, fitting these quan- 
tities does not give very strong constraints on the actual parameter 
values, but we can state the confidence level of any trends. For each 
of e, J, yin, and Lin, the deviation from NT as \h/r\ — » is less than 
5% with a confidence of 95%. For j integrated over ±2\h/r\, D[j] 
decreases with \h/r\ and increases with a/M both with 99% con- 
fidence. When integrating j over all angles, D[j] only decreases 
with \h/r\ to 99% confidence. For j^y, integrated over ±2\h/r\, D[ji„] 
only increases with a/M with 99.8% confidence and only decreases 
with \h/r\ with 97% confidence. When integrating over all an- 
gles, Z)[7in] only increases with a/M to essentially 100% confidence 
and only decreases with \h/r\ to 99.8% confidence. For e integrated 
over ±2\h/r\, D[e] only increases with \h/r\ with 98% confidence 
with no significant trend with a/M. When integrating e over all an- 
gles, D[e] only increases with a/M with 95% confidence with no 
significant trend with \h/r\. For Lin, there is a 98% confidence for 
this to increase with \h/r\ with no significant trend with a/M. Over- 
all, the most certain statement that can be made is that our results 
are strongly consistent with all deviations from NT becoming less 
than a few percent in the limit that \h/r\ —> across the full range 
of black hole spins. 



8 THIN DISKS WITH VARYING MAGNETIC FIELD 
GEOMETRY 

We now consider the effects of varying the initial field geometry. 
Since the magnetic field can develop large-scale structures that do 
not act like a local scalar viscosity, there could in principle be long- 
lasting effects on the accretion flow properties as a result of the 
initial field geometry. This is especially a concern for geometri- 
cally thin disks, where the 1-loop field geometry corresponds to 
a severely squashed and highly organized field loop bundle with 
long-range coherence in the radial direction, whereas our fiducial 
4-loop model corresponds to nearly circular loops which impose 
much less radial order on the MRI-driven turbulence. To investi- 
gate this question we have simulated a model similar to our fiducial 
run except that we initialized the gas torus with a 1-loop type field 
geometry instead of our usual multi-loop geometry. 

Figure |2T] shows the radial dependence of j, e, and T for 



the two field geometries under consideration, and Table |4] reports 
numerical estimates of various quantities at the horizon. Consider 
first the solid lines (4-loop fiducial run) and dotted lines (1-loop 
run) in Figure[2n both of which correspond to integrations in 9 over 
±2\h/r\ around the midplane. The simulation with 4-loops is clearly 
more consistent with NT than the 1-loop simulation. The value of j 
at the horizon in the 4-loop model deviates from NT by -2.9%. Be- 
tween the times of 12900M and 17300M, the 1-loop model deviates 
by -5.6%, while at late time over the saturated period the 1-loop 
model deviates by -7.2%. The long-dashed lines show the effect 
of integrating over all 6 for the 1-loop model. This introduces yet 
another systematic deviation from NT (as already noted in 35.7b : 
now the net deviation of j becomes -10.7% for times 12900M to 
17300M and becomes -15.8% for the saturated state. Overall, this 
implies that the assumed initial field geometry has a considerable 
impact on the specific angular momentum profile and the stress in- 
side the ISCO. This also indicates that the saturated state is only 
reached after approximately 17000M, and it is possible that the 1- 
loop model may never properly converge because magnetic flux of 
the same sign (how much flux is initially available is arbitrary due 
to the arbitrary position of the initial gas pressure maximum) may 
continue to accrete onto the black h ole and lead to a qual i tatively 
different accretion state (as seen in llgumenshchev et al.l ( l2003al) 
and[McKinnev & Gammie (200^ for their vertical field model). 
At early times, the nominal efficiency, e, shows no significant de- 
pendence on the field geometry, and sits near the NT value for both 
models. At late time in the 1-loop model, e rises somewhat, which 
may indicate the start of the formation of a qualitatively different 
accretion regime. 

Figure|22]shows the normalized luminosity. We see that the 1- 
loop model produces more luminosity inside the ISCO. For times 
12900M to 17300M, Li„ = 5.4% (integrated over all 9) compared 
to 3.5% for the 4-loop field geometry. Thus there is 50% more ra- 
diation coming from inside the ISCO in this model. At late time 
during the saturated state, Lin = 4.6% (integrated over all 9). Thus 
there is approximately 31% more radiation coming from inside the 
ISCO in this model during the late phase of accretion. 

Table|4]also reports the results for thick {\h/r\ x 0.3) disk mod- 
els initialized with the multi-loop and 1-loop field geometries. This 
again shows that the deviations from NT are influenced by the ini- 
tial magnetic field geometry and scale with \h/r\ in a way expected 
by our scaling laws. The I -loop models show deviations from NT 
in ] are larger as related to the larger value of T. The deviations 
from NT are less affected by the initial magnetic field geometry for 
thicker disks, because the deviations from NT are also driven by 
thermal effects and Reynolds stresses rather than primarily electro- 
magnetic stresses as for thin disks. 

These effects can be partially understood by looking at the 
specific electromagnetic stress, Y, shown in Figure |2T] We find 
T X 0.45 for the 4-loop field geometry. For times 12900M to 
17300M, Y X 0.71 in the 1-loop field geometry, and during the 
saturated state Y « 1.28. For times 12900M to 17300M, the 50% 
larger Y appears to be the reason for the 50% extra luminosity 
inside th e ISCO in the I-loop model. The magnetized thin disk 
model o f|G ammid ( 1 1 9991) predicts that, for a/M = 0, specific mag- 
netic fluxes of Y = 0.45, 0.71 should give deviations from NT of 
~Dlj] = 1.9, 3.9, respectively. These are close to the deviations 
seen in the simulations, but they are not a perfect match for rea- 
sons we can explain. First, the details of how one spatially-averages 
quantities (e.g., average of ratio vs. ratio of averages) when com- 
puting Y lead to moderate changes in its value, and, for integra- 
tions outside the midplane, comparisons to the Gammie model can 
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Table 4. Field Geometry Study 



Model Name \li/r\ M e D[e] j D[j] ji„ D[j,n] s Li„ Lei, lOOO,- T 
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Figure 21. Radial profiles of j and ji„ (top panel), e (middle panel), and 
T (bottom panel) are shown for two different initial field geometries. Re- 
sults for the fiducial 4-loop field geometry (model AOHR07) integrated over 
±2\h/r\ around the midplane ai'e shown by solid lines, for the 1-loop field 
geometry (model A0HR07LOOP1) integrated over ±2\h/r\ around the mid- 
plane by dotted lines, and the 1-loop model integrated over all angles by 
long-dashed lines. The short-dashed lines in the top two panels show the 
NT result. We see that the 1-loop field geometry shows larger deviations 
from NT in j and T compai'ed to the 4-loop geometry. The panels also 
reemphasize the point that including all 6 angles in the angular integration 
leads to considerable changes in j and T due to the presence of magnetic 
field in the corona-wind-jet. 



require slightly higher T than our diagnostic reports. Second, the 
finite thermal pressure at the ISCO leads to (on average over time) 
a deviation already at the ISCO that is non-negligible compared 
to the deviation introduced by electromagnetic stresses between 
the ISCO and horizon. This thermal component is not always im- 
portant, e.g., see the compar ison in Figure [TT] Still, as found in 
iMcKinnev & Gammid liooi) for thick disks at least, the deviations 
from NT contributed by the thermal pressure are of the same order 
as the deviations predi cted by the Gam mie model. These results 
motivate extending the iGammiel i ll 9991) model to include a finite 



Figure 22. Similar to Figure |2T] for the initial 4-loop and 1-loop field ge- 
ometries, but here we show the luminosity (top panel) and log-derivative 
of the luminosity (bottom panel). The luminosity is slightly higher for the 
1-loop model compared to the 4-loop model. 



(but still small) thermal pressure such that the boundary conditions 
at the ISCO lead to a non-zero radial velocity. 

Within the ISCO, we find that the time-averaged and volume- 
averaged comoving field strength for the 4-loop geometry roughly 



follows 1^1 



within ±2\h/r\ of the disk midplane, while 



at higher latitudes we have a slightly steeper scaling. For times 
I2900M to I7300M, the I-loop geometry has \b\ oc r"'-' within 
±2\h/r\ of the disk midplane, and again a slightly steeper scal- 
ing in the corona. Other than this scaling, there are no qualitative 
differences in the distribution of any c omoving field co mponent 
with height above the disk. While the iGammiel ( 1 19991) solution 
does not predict a power-law dependence for \b\, for a range be- 
tween T = 0.4-0.8, the variation near the horizon is approximately 
\b\ oc r""-^ - r""', which is roughly consistent with the simulation 
results. The slightly steeper slope we obtain for the 1-loop geome- 
try is consistent with a higher specific magnetic flux, although the 
variations in T for integration over different ran ges of angle im- 
ply stratification and a non-radial flow which the iGammid l fl999l) 
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Figure 23. Similar to Figure|7] but here we compare the initial 4-loop fidu- 
cial model (black solid lines) and the 1-loop model (dashed magenta lines). 
The horizontal black dashed lines in the second and third panels show the 
predictions of the NT model. The mass accretion rate, M, has larger root- 
mean-squared fluctuations in the 1-loop model, which is indicative of more 
vigorous turbulence. The nominal efficiency, e, shows no clear difference. 
The specific angular momentum, j, is lower in the 1-loop model compared 
to the 4-loop model. This indicates that the 1-loop field leads to larger stress 
within the ISCO. The absolute magnetic flux (per unit initial total absolute 
flux) on the black hole is larger in the 1-loop model than the 4-loop model. 
Since ~ 1/2 for the 1-loop model, essentially half of the initial loop 
was advected onto the black hole, while the other half gained angular mo- 
mentum and has been advected away. This may indicate that the 1-loop 
geometry is a poor choice for the initial field geometry, since the magnetic 
flux that ends up on the black hole is determined by the initial conditions. 
For times 12900M to 17300M, the value of T is about twice higher in the 
1-loop model, which implies about twice greater electromagnetic stresses 
within the ISCO. 



model does not account for. This fact and the rise in T with de- 
creasing radius seen in Figure [21] indicate a non-trivial degree of 
angular compression as the flow moves towards the horizon. Over- 
all, our results suggest that deviati ons from NT de pend on the as- 
sumed field geometry, and that the lGamrni3 ( 1 19991) model roughly 
fits the simulations. 

Figure [23] shows the same type of plot as in Figure [7] but 
here we compare the fiducial 4-loop model with the 1-loop model. 
As mentioned above, the 1-loop geometry has a larger deviation 
in } from the NT value, corresponding to larger stresses inside 
the ISCO. The absolute magnetic flux (per unit initial total ab- 
solute magnetic flux) on the black hole <l>r is of order 1 /2, sug- 
gesting that the inner half of the initial field bundle accreted onto 
the black hole, while the other half was advected to larger radii. 
This is consistent with what i s seen in simulations of thick tori 
jMcKinnev & Gammid |2004| : iBeckwith et al] l2008al) . This sug- 
gests that using the 1-loop geometry leads to results that are sen- 
sitive to the initial absolute magnetic flux, while the multiple loop 
geometry leads to results that are insensitive to the initial abso- 
lute magnetic flux. Such dependence of the electromagnetic stress 
on initial magnetic field geometry has also been reported in 3D 



Figure 24. Normalized electromagnetic stress, W/M, as a function of ra- 
dius for the fiducial model (black lines) and the otherwise identical 1-loop 
model (magenta lines). The solid lines correspond to a 6 integration over 
all angles, while the dotted lines correspond to a S integration over ±2|/i/r|. 
The dashed fine shows the NT result, while the dashed green lines show 
the Gammie (1999) result for T = 0.60,0.89,0.90, 1.21 for lines from the 
bottom to the top. The Gammie (1999) model gives a reasonable fit to the 
simulation's stress profile within the ISCO. The 1-loop model shows about 
50% larger peak normalized stress (integrated over all angles) compared to 
the multi-loop model (integrated over all angles), which is consistent with 
the 1-loop model leading to larger deviations from NT (about 50% larger 
luminosity over the period used for time averaging). The large differences 
between the solid and dotted lines again highlights the fact that the stress 
within the disk is much smaller than the stress over all 8 that includes the 
corona+wind+jet. As pointed out in SOS, even though such a plot of the 
electromagnetic stress appears to indicate large deviations from NT within 
the ISCO, this is misleading because one has not specified the quantitative 
eft'ect of the non-zero value of WjM on physical quantities within the ISCO. 
Apparently large values of Wj M do not necessarily correspond to large de- 
viations from NT. For example, quantities such as j, e. and the luminosity 
only deviate by a few percent from NT for the multi-loop model. 



pseudo-Newtonian simulations by Revnolds & Armitagel ( 1200 ll) 
and in 3D GRMHD simulations bv lBeckwith et al.l ( l2008a[) . 

Figure l24l shows the electromagnetic stress as computed by 
equation J27t for the multiple loop fiducial model (A0HR07) and 
the otherwise identical 1-loop model (A0HR07LOOP1). We only 
show the electromagnetic part of the stress, and within the ISCO 
this is to within a few percent the same as the total stress ob- 
tained by including all terms in the stress-energy tensor. Out- 
side the ISCO, the total stress agrees more with the NT model. 
The figure shows the full-angle integrated electromagnetic stress, 
the electromag netic stre s s integ rated over only ±2\h/r\, the NT 
stress, and the iGammiei j 19991) electromagnetic stress for T = 
0.60, 0.89, 0.90, 1.21 (we choose T, the only free parameter of the 
model, such that the peak magnitude of the stress agrees with the 
simulation). The chosen T values are close to our d iagnostic's value 
of T for these models, which demonstrates that the lGamrni 3i ll999l) 
model is consistently predicting the simulation's results with a sin- 
gle free parameter. The stress is normalized by the radially depen- 
dent M(r) that is computed over the same 6 integration range. We 
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do not restrict the integration to bound material as done in SOS 
(in SOS, thie stress is integrated over ±2\h/r\ and only for bound 
material, while in NIO the stres|3 is only over bound material). 
The stress for the fiducial model was time-averaged over the sat- 
urated state, while the 1-loop model was time-averaged from time 
12900M to 17300M in order to best compare with the early phase 
of accretion for the 1-loop model studied in NIO. 

Figurel24lshows that the simulation and NT stress do not agree 
well, and it suggests there is an apparently large stress within the 
ISCO. However, as first pointed out by SOS and discussed in sec- 
tion [3]4l this stress does not actually correspond to a large devia- 
tion from NT in physically relevant quantities such as the specific 
angular momentum, specific energy, and luminosi t y. Thi s point is 
clarified by making a comparison to the iGanmii model's 
stress, which agrees reasonably well with the simulation stress in- 
side the ISCO. Even though the stress may appear large inside the 
ISCO, the stress corresponding to the Gammie model with (e.g.) 
Y = 0.60 only translates into a few percent deviations from NT. 
This figure also demonstrates that the initial magnetic field geom- 
etry affects the amplitude of the stress in the same direction as 
it affects other quantities and is reasonably well predicted by the 
lGammiel ( fT99^ model. The initial magnetic field sets the saturated 
value of Y, which is directly related to the electromagnetic stresses 
within the ISCO. The 1-loop model leads to a peak stress (inte- 
grated over all angles) within the ISCO that is about 50% larger 
than the multi-loop model (integrated over all angles), which is 
likely related to the extra 50% luminosity in the 1-loop model com- 
pared to the multi-loop model. The fact that the stress normaliza- 
tion changes with initial field geometry is consistent w i th othe r 
3D GRMHD simulations of thick disks bv lSeckwith et alUlOOSah . 
This figure again shows how the stress within the disk (±2\h/r\) is 
much smaller than the total disk+corona+wind+jet (all 9). 

Finally, we discuss previous results obtained for other field 
geom etries using an energy-conserv ing two-dimensional GRMHD 
code jMcKinnev & Gammie|[2004h . While such two-dimensional 
simulations are unable to sustain turbulence, the period over which 
the simulations do show turbulence agrees quite well with the cor- 
responding period in three-dimensional simulations. This implies 
that the turbulent period in the two-dimensio nal simulations may be 
qualita tively correct. The fiducial model of iMcKinnev & Gammid 
( |2004|) was of a thick (\h/r\ ~ 0.2-0.25) disk with a 1-loop initial 
field geometry around an ajM = 0.9375 black hole. This model 
had Y ~ 1 near the midplane within the ISCO and Y ~ 2 when 
integrated over all Q angles. Their measured value of j « 1.46 in- 
tegrated o ver all angles, \b\ oc within the ISCO within the disk 
midplane jMcRinnev & N aravan 2007b), along with Y ~ 1-2 are 
roughly consistent with the Gammie (199^ model prediction of 
] X 1.5. Similarly, the strong vertical field geometry model they 
studied had Y ~ 2 near the midplane within the ISCO and Y ~ 6 in- 
tegrated over all 6 angles. Their measurement of y ~ - 1 integrated 
over all angles is again roughly consistent with the model predic- 
tion of y x -1.2 for Y ~ 6. Note that in this model, Y rises (as usual 
to reach saturation) with time, but soon after Y > 2 in the midplane, 
the disk is pushed away by the black hole and then Y is forced to 
be even larger. Evidently, the accumulated magnetic flux near the 
black hole pushes the system into a force-free magnetosphere state 



NlO's figures 12 and 13 show stress vs. radius, but some of the integrals 
they computed were not re-normalized to the full 2n when using their sim- 
ulation i^-box size of njl, so their stress curves are all a constant factor of 4 
times larger than the actual stress (Noble, private communication). 



- not an accretion state. This shows the potential danger of using 
strong-field initial conditions (like the 1-loop field geometry), since 
the results are sensitive to the assumed initial flux that is placed 
on (or rapidly drops onto) the black hole. Even while the disk is 
present, this particular model exhibits net angular momentum ex- 
traction from the black hole. This interesting result needs to be con- 
firmed using three-dimensional simulations of both thick and thin 
disks. 



9 COMPARISONS WITH OTHER RESULTS 

The resul ts we have ob tained in the present work are consistent with 
those of lArmitageTt al. (2001) and Reynolds & Fabian (200^, 
who carried out pseudo-Newtonian studies, and with the results 
of SOS, who did a full GRMHD simulation. Both of these studies 
found only minor deviations from NT for thin accretion disks with 
a multi-loop initial field geometry. However, more recently, N09 
and NIO report apparently inconsistent results, including factors of 
up to five larger deviations from NT in the specific angular momen- 
tum (2% in SOS versus 10% in NIO) for the same disk thickness of 
\hlr\ ~ 0.07. They also find a 50% larger deviation from NT in 
the luminosity (4% in SOS versus 6% in N09). Furthermore, in NIO 
they conclude that the electromagnetic stresses have no dependence 
on disk thickness or initial magnetic field geometry, whereas we 
find that the electromagnetic stresses have a statistically significant 
dependence on both disk thickness and magnetic field geometry. 

We have considered several possible explanations for these 
differences, as we now describe. We attempt to be exhaustive in 
our comparison with the setup and results by N09 and NIO, be- 
cause our works and their works seek accuracies much better than 
order two in measuring deviations from NT. Thus, any deviations 
between our results by factors of two or more must be investigated 
further in order to ensure a properly understood and accurate result. 

First, we briefly mention some explanations that NIO propose 
as possible reasons for the discrepant results, viz., differences in I) 
numerical algorithm or resolution; 2) box size in ^-direction: Ai^; 
3) amplitude of initial perturbations; 4) accuracy of inflow equilib- 
rium; and 5) duration of the simulations. Our algorithms are sim- 
ilar except that their PPM interpolation scheme assumes primitive 
quantities are cell averages (standard PPM), while ours assumes 
they are point values (as required to be applied in a higher-order 
scheme). They used LAXF dissipative fluxes, while we used HLL 
fluxes that are about twice more accurate for shocks and may be 
more accurate in general. On the other hand, they used parabolic 
interpolation for the Toth electric field, while we use the standard 
Toth scheme. Given these facts, we expect that the accuracy of our 
algorithms is similar. Overall, our convergence testing and other di- 
agnostics (see ^ confirm that none of their proposed issues can be 
the cause of differences between SOS and NIO. 

We have shown that inflow equilibrium must include satu- 
ration of the specific magnetic flux, Y, which generally saturates 
later in time than other quantities. By running our fiducial model 
A0HR07 to a time of nearly 30000M, we ensure that we have a 
long period of steady state conditions to compute our diagnostic 
quantities. The fact that we need to run our fiducial thin disk sim- 
ulation for such a long time to reach inflow equilibrium up to a ra- 
dius r ~ 9M is completely consistent with our analytical estimate 
of the time scale calculated using Eq. IB7I of Appendix |B] (see the 
earlier discussion in 33.5| and Fig.|6j. In the comparison between 
the numerical and analytical results shown in Figure (6] we found 
agreement by setting a\h/r\^ a 0.00033 which, for our disk with 
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\hlr\ a; 0.064, corresponds to cf a 0.08. With this value of a\h/r\^, 
we would have to run the simulation until t ~ 83000M, 160000M, 
to reach inflow equilibrium out to 15M, 20M, respectively, corre- 
sponding to a couple viscous timescales at that radius. NIO state 
that they reach inflow equilibrium within r ~ 15M-20M in a time 
of only / ~ lOOOOM. Since their disk thickness is \h/r\ a; 0.06, even 
a single viscous timescale would require their simulations to have 
a ~ 0.38 to reach inflow equilibrium up to r ~ 15M, and an even 
larger value of o- for r ~ 20M. This seems unlikely. We can partially 
account for their result by considering our 1-loop model, which up 
to t ~ 17000M has a\h/r\- twice as large and a about 70% larger 
than in the fiducial 4-loop run. However, this still falls far short by 
a factor of roughly 3 of what NIO would require for inflow equilib- 
rium up to 15M-20M. Further, our A0HR07LOOP1 model, which 
is similar to their model, only reaches a saturated state by 17000M, 
and only the Y quantity indicates that saturation has been reached. 
If we were to measure quantities from lOOOOM to 15000M as in 
NIO, we would have underestimated the importance of magnetic 
field geometry on the electromagnetic stresses. 

Since all these simulations are attempting to obtain accuracies 
better than factors of two in the results, this inflow equilibrium issue 
should be explored further. A few possible resolutions include: 1) 
NlO's higher resolution leads to a much larger a; 2) their disk has a 
larger "effective" thickness, e.g. \h/r\ ~ 0.13, according to Eq. 5.9.8 
in NT (see Eq. lB7l of Appendix|B}; 3) some aspects of their solution 
have not yet reached inflow equilibrium within a radius much less 
than r ~ 15M, such as the value of Y vs. time that saturates much 
later than other quantities; or 4) they achieve constant fluxes vs. 
radius due to transient non- viscous effects - although one should be 
concerned that the actual value of such fluxes continues to secularly 
evolve in time and one still requires evolution on the longer viscous 
(turbulent diffusion) timescale to reach true inflow equilibrium. 

Second, we considered various physical setup issues, includ- 
ing differences in: 1) range of black hole spins considered; 2) range 
of disk thicknesses studied; 3) ad hoc cooling function; and 4) equa- 
tion of state. We span the range of black hole spins and disk thick- 
nesses studied by NIO, so this is unlikely to explain any of the 
differences. Some differences could be due to the disk thickness 
vs. radius profiles established by the ad hoc cooling functions in 
the two studies. NlO's cooling function is temperature-based and 
allows cooling even in the absence of any dissipation, while ours 
is based upon the specific entropy and cools the gas only when 
there is dissipation. Both models avoid cooling unbound gas. In 
S08 and in the present paper, we use an ideal gas equation of state 
with r = 4/3, while N09 and NIO used T = 5/3. The proper- 
ties of the turbulence do appe ar to depend on the equation of state 
jMignone & McKinnevll2007h . so it is important to investigate fur- 
ther the role of T in thin disks. 

Third, the assumed initial field geometry may introduce crit- 
ical differences in the results. Issues with the initial field geome- 
try include how many field reversals are present, how isotropic the 
field loops are in the initial disk, how the electromagnetic field en- 
ergy is distributed in the initial disk, and how the magnetic field is 
normalized. In SOS and here, we have used a multi-loop geometry 
in the initial torus consisting of alternating polarity poloidal field 
bundles stacked radially. We ensure that the field loops are roughly 
isotropic within the initial torus. We set the ratio of maximum gas 
pressure to maximum magnetic pressure to /?maxcs = 100, which 
gives us a volume-averaged mean f} within the dense part of the 
torus (po/po,mdx > 0.2) oip~ 800. Our procedure ensures that all 
initial local values of f} within the disk are much larger than the val- 



ues in the evolved disk, i.e., there is plenty of room for the magnetic 
field to be amplified by the MRI. 

We have also studied a 1-loop geometry that is similar to the 
1-loop geometry used in N09 and NIO. Their initial (/i-component 
of the vector potential is oc MAX(po/po,max - 0.25,0) (Noble, 
private communication). They initialize the magnetic field geome- 
try by ensuring that the volume-averaged gas pressure divided by 
the volume-averaged magnetic pressure is /3ai,crages = 100 (Noble, 
private communication). (They stated that the mean initial plasma 
/3 is yS = 100.) For their thin disk torus model parameters, this nor- 
malization procedure leads to a portion of the inner radial region of 
the torus to have a local value of /? ~ 3 - 8, which may be a source 
of differences in our results. Such a small p is lower than present in 
the saturated disk. NIO make use of an old er set of simulations from 
a different non-en ergy-conserving code l lHawlev&Krolikll2006l ; 
iBeckwith et al HoOSaD to investigate the effect of other field geome- 
tries. The results from this other code have strong outliers, e.g., the 
KDOc model, and so we are unsure if these other simulations can 
be used for such a study. 

NIO state that they find no clear differences in the electro- 
magnetic stresses for different initial field geometrie s. As shown 
in their figures 12 and 13, the I Agol & KrolikI bOOd) model cap- 
tures the smoothing of the stress outside the ISCO, but it is not 
a model for the behavior of the stress inside the ISCO. We find 
that electromagnetic stresses have a clear dependence on both disk 
thickness and t he initial magnet ic field geometry, with a trend that 
agrees with the iGammi model of a magnetized thin disk. 

Our Figure[24lshows that the stress w ithin the ISCO is reasonably 
well modelled by the lGammid j 1 9991) model. Our 1-loop thin disk 
model gives a peak normalized stress (integrated over all angles) of 
about 3.2x 10"^ for times 12900M to 17300M, which is comparable 
to the I -loop thin disk model in NIO with peak normalized stress 
(integrated over all angles) of about 2.5 x 10"^ (after correcting for 
their (f>-box size). Hence, we are able to account for the results of 
their 1-loop model. 

In addition, we used the specific magnetic flux, Y, an ideal 
MHD invariant that is conserved within the ISCO, to identify how 
electromagnetic stresses scale with disk thickness and magnetic 
field geometry. In the saturated state, the value of Y, which con- 
trols the electromagnetic stresses, is different for different initial 
magnetic field geometries. We find that j within the disk (±2|/i/r| 
from the midplane) deviates from NT by -3% in our 4-loop model 
and -6% in our I-loop model for times 12900M to 17300M. Inte- 
grating over all angles, j deviates by -6% for the 4-loop model and 
-11% for the 1-loop model for times 12900M to 17300M. Thus, 
we find a clear factor of two change, depending on the assumed 
initial field geometry and the range of integration. The excess lu- 
minosity is 3.5% for the 4-loop model and 5.4% for the 1-loop 
model for times 12900M to 17300M. Recalling that NIO find a 
deviation from NT of about -10% in j (integrated over all angles) 
and a luminosity excess beyond NT of about 6%, this shows we can 
completely account for the apparent inconsistencies mentioned by 
NIO by invoking dependence of the results on the initial field ge- 
ometry and the presence of extra stress beyond the disk component 
of the accretion flow. 

Fourth, let us consider measurement and interpretation differ- 
ences. Our ultimate goal is to test how well the NT model describes 
a magnetized thin accretion disk. The primary quantity that is used 
to measure this effect in S08 and NIO is the specific angular mo- 
mentum J. However, the measurements are done differently in the 
two studies. In S08 as well as in this paper, we focus on the disk 
gas by limiting the range of 6 over which we compute the averag- 
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ing integrals {±2\h/r\ from the midplane). In contrast, NIO compute 
their integrals over a much wider range of 9 which includes both 
the disk and the corona-wind-jet (Noble, 2010, private communi- 
cations). We have shown in § l5.7l that the disk and corona-wind-jet 
contribute roughly equally to deviations of j from the NT value. 
In principle, the luminosity from the corona-wind-jet could be im- 
portant, but we have shown that the excess luminosity of bound gas 
within the ISCO is dominated by the disk. This means that the mea- 
sure used by NIO, consisting of integrating over all gas to obtain j, 
cannot be used to infer the excess luminosity of bound gas within 
the ISCO. Further, the corona would largely emit non-thermal radi- 
ation, so for applications in which one is primarily interested in the 
thermal component of the emitted radiation, one should evaluate 
the accuracy of the NT model by restricting the angular integration 
range to the disk component within ±2\h/r\. 

Fifth, let us consider how the results from NIO scale with 
disk thickness for the specific case of a non-spinning (a/M = 0) 
black hole. We have performed a linear least squares fit of their 
simulation results, omitting model KDOc which is a strong out- 
lier. For / integrated over all 6, their relative difference follows 

a; -7 - 45|/i/r| with confidence of 95% that these coefficients, 
respectively, only deviate by ±67% and ±89%. These fits imply 
that, as \h/r\ — > 0, the relative deviation of j from the NT value is 
about -7%, but they could easily be as low as -2%. Their results 
do not indicate a statistically significant large deviation from NT 
as \h/r\ — > 0. Since the total deviation in j from NT includes the 
effects of electromagnetic (and all other) stresses, this implies that 
their models are consistent with weak electromagnetic stresses as 
\h/r\ 0. 

Further, we have already established that the 1-loop geometry 
gives (at least) twice the deviation from NT compared to the 4-loop 
geometry, plus there is another factor of two arising from including 
the corona- wind-jet versus not including it. This net factor of 4 
applied to NlO's results implies that j would deviate by about -2% 
or even as low as -0.5% from NT in the limit \h/r\ — > if they were 
to consider a 4-loop field geometry and focus only on the disk gas. 
Thus, their models show no statistically significant large deviations 
from NT. In addition, our results in section ItTI show that, whether 
we consider an integral over all angles or only over the disk, there 
is no statistically significant large deviation from NT as \h/r\ — > 0. 

In summary, we conclude that the apparent differences be- 
tween the results obtained in SOS and the present paper on the one 
hand, and those reported in N09 and NIO on the other, are due to 
I) dependence on initial magnetic field geometry (multi-loop vs 
I -loop); 2) dependence upon the initial magnetic field distribution 
and normalization; and 3) measurement and interpretation differ- 
ences (disk vs. corona- wind-jet). Note in particular that the 1-loop 
initial field geometry is severely squashed in the vertical direction 
and elongated in the radial direction for thin disks, and it is not clear 
that such a geometry would ever arise naturally. There are also indi- 
cations from our simulation that the I-loop geometry may actually 
never reach a converged state due to the arbitrary amount of mag- 
netic flux accreted onto the black hole due to the single polarity of 
the initial magnetic field. Finally, if one is trying to test how well 
simulated thin accretion disks compare with NT, then it is impor- 
tant to restrict the comparison to disk gas near the midplane. One 
should not expect the gas in the corona-wind-jet to agree with the 
NT model. 



10 DISCUSSION 

We now discuss some important consequences of our results and 
also consider issues to be addressed by future calculations. First, 
we discuss the relevance to black hole spin measurements. 

In recent years, black hole spin parameters have been mea- 
sured in several black hole x-ray binaries by fitting their x-ray 
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method is based o n several assumptions that require testing 
jNaravan et alj|2008h . the most critical being the assumption that 
an accretion disk in the radiatively-efficient thermal state is well- 
described by the Novikov-Thome model of a thin disk. More 
specifically, in analyzing and fitting the spectral data, it is assumed 
that the radial profile of the radiative flux, or equivalently the effec- 
tive temperature, in the accretion disk closely follows the prediction 
of the NT model. 

Practitioners of the continuum-fitting method generally re- 
strict their attention to relatively low-luminosity systems below 
30% of the Eddington luminosity. At these luminosities, the max- 
imum height of the disk photosphere above th e midplane is less 
than 10% of the radius, i.e., (h / r)photoiphas: ^ O.I jMcClintock et al.l 
l2006l) . For a typical disk, the photospheric disk thickness is ap- 
proximately twice the mean absolute thickness \h/r\ that we con- 
sider in this paper. Therefore, the disks that observers consider for 
spin measurement have \h/r\ < 0.05, i.e., they are thinner than the 
thinnest disk (\h/r\^i„ ~ 0.06) that we (SOS, this paper) and others 
(N09, NIO) have simulated. 

The critical question then is the following: Do the flux pro- 
files of very thin disks match the NT prediction? At large radii 
the two will obviously match very well since the flux profile is 
determined simply by energy conservatiorF^. However, in the re- 
gion near and inside the ISCO, analytic models have to apply a 
boundary condition, and the calculated flux profile in the inner re- 
gion of the disk depends on this choice. The conventional choice 
is a "zero-torque" boundary condition at the ISCO. Unfortunately, 
there is disagreement on the validity of this assumption. Some 
authors have argued that the magnetic field strongly modifies the 
zero-torque condition and that, therefore, real disks might behave 
very differentl y from the NT model near the ISCO (iKrolik I999I; 
lGammielll999l) . Other authors, based either on heuristic arguments 
or on hydrodynamic calculations, find that the NT model is accu- 
rate even near th e ISCO so long as the disk i s geometrically thin 



(Paczvnski 200d ;lAfshordi & Paczvnskill2003l ; IShafee et al.ll200S j ; 
lAbramowicz et aij|201(ll) . Investigating this question was the pri- 
mary motivation behind the present study. 

We described in this paper GRMHD simulations of geometri- 
cally thin (\h/r\ ~ 0.07) accretion disks around black holes with a 
range of spins: a/M = 0, 0.7, 0.9, 0.98. In all cases, we find that 
the specific angular momentum j of the accreted gas as measured 
at the horizon (this quantity provides information on the dynamical 
importance of torques at the ISCO) shows only minor deviations 
at the level of ~ 2%^% from the NT model. Similarly, the lumi- 
nosity emitted inside the ISCO is only ~ 3%-7% of the total disk 
luminosity. When we allow for the fact that a large fraction of this 
radiation will probably be lost into the black hole because of rel- 
ativistic beaming as the gas plunges inward (an effect ignored in 

This is why the formula for the flux as a function of radius in the standard 
thin disk model d oes not depend on details like the viscosity parameter a 
jpranket all 19921) . 
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our luminosity estimates), we conclude that the region inside the 
ISCO is likely to be quite unimportant. Furthermore, our investi- 
gations indicate that deviations from the NT model decrease with 
decreasing \hlr\. Therefore, since the disks of interest to observers 
are generally thinner than the thinnest disks we have simulated, the 
NT model appears to be an excellent approximation for modeling 
the spectra of black hole disks in the thermal state. 

One caveat needs to be mentioned. Whether or not the total 
luminosity of the disk agrees with the NT model is not important 
since, in spectral modeling of data, one invariably fits a normal- 
izatio n (e.g., the accretion rate M in the model KERRBB: lLi et all 
which absorbs any deviations in this quantity. What is impor- 
tant is the shape of the flux profile versus radius. In particular, one 
is interested i n the radius at which the flux or effective tem perature 
IS maximum jNaravan et alJl2008l ; lMc"ciintock et al1l2009l) . Quali- 
tatively, one imagines that the fractional shift in this radius will be 
of order the fractional torque at the ISCO, which is likely to be of 
order the fractional error in j. We thus guess that, in the systems of 
interest, the shift is nearly always below 10%. We plan to explore 
this question quantitatively in a future study. 

Another issue is the role of the initial magnetic field topol- 
ogy. We find that, for aJM = 0, starting with a 1-loop field ge- 
ometry gives an absolute relative deviation in / of 7.1%, and an 
excess luminosity inside the ISCO of 4.9%, compared to 2.9% and 
3.5% for our standard 4-loop geometry. Thus, having a magnetic 
field distribution with long-range correlation in the radial direction 
seems to increase deviations from the NT model, though even the 
larger effects we find in this case are probably not a serious con- 
cern for black hole spin measurement. Two comments are in order 
on this issue. First, the 4-loop geometry is more consistent with 
nearly isotropic turbulence in the poloidal plane and, therefore, in 
our view a more natural initial condition. Second, the I-loop model 
develops a stronger field inside the ISCO and around the black hole 
and might therefore be expected to produce a relativistic jet with 
measurable radio emission. However, it is well-known that black 
hole x-ray binaries in the thermal state have no detectable radio 
emission. This suggests that the magnetic field is probably weak, 
i.e., more consistent with our 4-loop geometry. 

Next, we discuss the role of electromagnetic stresses on the 
dynamics of the gas in the plunging region inside the ISCO. In 
order to better understand this issue, we have extracted for each 
of our simulations the radial profile of the specific magnetic flux, 
Y. This quantity appears as a dimensionless free parameter (called 
Fg ^) in the s i mple MHD model of the plunging region developed 
by lGamrni3 ( Il999l) . The virtues of the specific magnetic flux are 
its well-defin ed normalization and it s constancy with radius for sta- 
tionary flows l lTakahashi et aljl99(lh . In contrast, quantities like the 
stress W or the viscosity parameter a have no well-defined normal- 
ization; W can be normalized by any quantity that has an energy 
scale, such as po, M, or b^, while a could be defined with respect to 
the total pressure, the gas pressure, or the magnetic pressure. The 
numerical values of IV or or inside the ISCO can thus vary widely, 
depending on which definition one chooses. For instance, although 
SOS found a ~ \ inside the ISCO, the specific angular momentum 
flux, ], deviated from NT by no more than a few percent. Further, 
Figure [24lshows that (even for the multi-loop model) the stress ap- 
pears quite large within the ISCO, but this is misleading because 
the effects of the stress are manifested in the specific angular mo- 
mentum, specific energy, and luminosity - all of which agree with 
NT to within less than 10% for the multi-loop model. Since W and 
a do not have a single value within the ISCO or a unique normal- 



ization, we conclude that they are not useful for readily quantifying 
the effects of the electromagnetic stresses within the ISCO. 

Gammie's (1999) model shows how the value of T is directly 
related to the electromagnetic stresses within the ISCO. Unfortu- 
nately, the actual value of Y is a free parameter which cannot be 
easily determined from first principles. It is possible that accretion 
disks might have Y » I, in which case, the model predicts large de- 
viations from NT. For example, if T = 6, then for an ajM = black 
hole ] is lowered by 56% relative to the NT model. We have used 
our 3D GRMHD simulations which include self -consistent MRI- 
driven turbulence to determine the value of Y for various black hole 
spins, disk thicknesses, and field geometries. For the multiple-loop 
field geometry, we find that the specific magnetic flux varies with 
disk thickness and spin as 

h 



Y K 0.7 + 



-0.6i^, 



(41) 



within the disk component, which indicates that electromagnetic 

stresses are weak and cause less than 8% deviations in j in the limit 

— » for all black hole spins. Our rough analytical arguments 

for how Y should scale with \hlr\ and fl/M are consistent with the 

above formula. Even for the 1-loop field geometry, Y < 1 for thin 

disks, so electromagnetic stresses cause only minor deviations from 

NT for all black hole spins (for Y < 1, less than 12% in /). Not all 
I 1 J — ^ — I 

aspects of the Gammie (1999) model agree with our simulations. 
As found in McKinney & Gammie (2004), the nominal efficiency, 
e, does not match well and for thin disks is quite close to NT. Since 
the true radiative efficiency is limited to no more than e, this pre- 
dicts only weak deviations from NT in the total luminosity even 
if ] has non-negligible deviations from NT. Also, this highlights 
that the deviations from NT in j are due to non-dissipated electro- 
magnetic stresses and cannot be used to directly predict the excess 
luminosity within the ISCO. The assumption of a radial flow in a 
split-monopole field is approximately valid, but the simulations do 
show some non-radial flow and vertical stratification, a non-zero ra- 
dial velocity at the ISCO, and thermal energy densities comparable 
to magnetic energy densities. Inclusion of these effects is required 
for better consistency with simulation results inside the ISCO. 

Next, we consider how our results lend some insight into the 
spin evolution of black holes. Standard thin disk theory with pho- 
ton capture predicts that an accreting black hole spins up until it 
reaches spin equilibrium at a^^/M x 0.998 ( lThorne|[T974E . On the 
other hand, thick non-radiative accretion flows deviate significantly 
from NT and reach equilibrium at a^q/M ~ 0.8 for a model with 
a ~ 0.3 and \h/r\ ~ 0.4 near the horizon jPopham & Gamnii3 
Il998l) . GRMHD simulations of moderately thick {\h/r\ ~ 0.2- 
0.25) magnetized accretion flows give flcq/A^ ~ 0.9 jGammie et al.l 
|2004 . In this paper, we find from our multi-loop field geometry 
models that spin equilibrium scales as 

1.1 -0.8^ 

M 



(42) 



where one should set a^^/M = 1 if the above formula gives a^-^/M > 
I. This gives a result consistent with the above-mentioned studies 
of thick disks, and it is also consistent with our rough analytical pre- 
diction based upon our scaling of Y and using the Gammie model 
prediction for the spin equilibrium. This result also agrees with the 
NT result in the limit \h/r\ — > within our statistical errors, and 
shows that magnetized thin disks can approach the theoretical limit 
of a^q/M x 1, at least in the multi-loop case. In the single-loop field 
geometry, because of the presence of a more radially-elongated ini- 
tial poloidal field, we find a slightly stronger torque on the black 
hole. However, before a time of order 17000M, the deviations in 
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the equilibrium spin parameter, a^q/M, between the 4-loop and 1- 
loop field geometries appear to be minor, so during that time the 
scaling given above roughly holds. Of course, it is possible (even 
likely) that radically dilferent field geometries or anomalously large 
initial field strengths will lead to a different scaling. 

Lastly, we mention a number of issues which we have ne- 
glected but are potentially important. A tilt between the angular 
momentum vector of the disk and the blac k hole rotation axi s 
might significantly affect the accretion flow ^Fragile et alj|2007h . 
We have not accounted for any radiative transfer physics, nor 
have we attempted to trace photon trajectories (see, e.g. N09 and 
iNoble & Kro"iiMl2009l) . In principle one may require the simula- 
tion to be evolved for hundreds of orbital times at a given radius 
in ord er to completely erase the initial conditions jSorathia et al.l 
whereas we only run the model for a couple of viscous time 
scales. New pseudo-Newtonian simulations show that convergence 
may require resolving se veral disk scale heights with high resolu- 
tion jSorathia etai]|2010h . and a similar result has been found also 
for shearing box calculations with no net flux and no stratification 
(Stone 2010, private communication). In contrast, we resolve only 
a couple of scale heights. Also, we have only studied two differ- 
ent types of initial field geometries. Future studies should consider 
whether alternative field geometries change our results. 



11 CONCLUSIONS 

We set out in this study to test the standard model of thin 
accretion disks around rotating black holes as developed by 
iNovikov & Thomd ( Il973h . We studied a range of disk thicknesses 
and black hole spins and found that magnetized disks are consistent 
with NT to within a few percent when the disk thickness \h/r\ < 
0.07. In addition, we noted that deviations from the NT model de- 
crease as \h/r\ goes down. These results suggest that bl ack hole spin 
measurements via the x- ray continuum-fitting method llZhang et al.l 



I997':'Shafee et al.'2006VMcClintoc k et alj2006l ; lDavis et alj2006l ; 
Liu et al. 2008; Gou et al. 2009, 201^, which are based on the NT 



model, are robust to model uncertainties so long as \h/r\ < 0.07. At 
luminosities below 30% of Eddington, we estimate disk thicknesses 
to be \h/r\ ^ 0.05, so the NT model is perfectly adequate. 

These results were obtained by performing global 3D 
GRMHD simulations of accreting black holes with a variety of disk 
thicknesses, black hole spins, and initial magnetic field geometries 
in order to test how these affect the accretion disk structure, angular 
momentum transport, luminosity, and the saturated magnetic field. 
We explicitly tested the convergence of our numerical models by 
considering a range of resolutions, box sizes, and amplitude of ini- 
tial perturbations. As with all numerical studies, future calculations 
should continue to clarify what aspects of such simulations are con- 
verged by performing more parameter space studies and running 
the simulations at much higher resolutions. For example, it is pos- 
sible that models with different black hole spins require more or 
less resolution than the a = models, while we fixed the resolution 
for all models and only tested convergence for the a = models. 

We confirmed previous results by SOS for a non-spinning 
{a/M = 0) black hole, which showed that thin (\h/r\ < 0.07) disks 
initialized with multiple poloidal field loops agree well with the 
NT solution once they reach steady state. For the fiducial model 
described in the present paper, which has similar parameters as the 
SOS model, we find 2.9% relative deviation in the specific angular 
momentum accreted through the disk, and 3.5% excess luminosity 
from inside the ISCO. Across all black hole spins that we have con- 



sidered, viz., a/M = 0, 0.7, 0.9, 0.98, the relative deviation from 
NT in the specific angular momentum is less than 4.5%, and the 
luminosity from inside the ISCO is less than 7% (typically smaller, 
and much of it is likely lost to the hole). In addition, all deviations 
from NT appear to be roughly proportional to \h/r\. 

We found that the assumed initial field geometry modifies the 
accretion flow. We investigated this effect by considering two dif- 
ferent field geometries and quantified it by measuring the specific 
magnetic flux, T, which is an ideal MHD invariant (like the spe- 
cific angular momentum or specific energy). The specific magnetic 
flux can be written as a dimension less free parame ter that enters 
the magnetized thin disk model of iGamrnig jl999l) . This param- 
eter determines how much the flow deviates from NT as a result 
of electromagnetic stresses. We found that T allows a quantita- 
tive understanding of the flow within the ISCO, while the electro- 
magnetic stress (W) has no well-defined normalization and varies 
widely within the ISCO. While a plot of the stress may appear to 
show large stresses within the ISCO, the actual deviations from NT 
can be small. This demonstrates that simply plotting W is not a use- 
ful diagnostic for measuring deviations from NT. We found that the 
specific magnetic flux of the gas inside the ISCO was substantially 
larger when we used a single poloidal magnetic loop to initialize the 
simulation compared to our fiducial 4-loop run. For a/M = and 
\h/r\ < 0.07, the eariy saturated phase (times 12900M to 17300M) 
of the evolution for the I-loop geometry gave 5.6% relative de- 
viation in the specific angular momentum and 5.8% excess lumi- 
nosity inside the ISCO. These deviations are approximately twice 
as large as the ones we found for the 4-loop simulation. At late 
times, the 1-loop model generates significant deviations from NT, 
which is a result similar to th at found in a vertical field model in 
iMcKinnev & Gammid hoo4) . However, we argued that the mul- 
tiple loop geometry we used is more natural than the single loop 
geometry, since for a geometrically thin disk the magnetic field in 
the 1-loop model is severely squashed vertically and highly elon- 
gated radially. The 1-loop model is also likely to produce a strong 
radio jet. 

More significant deviations from NT probably occur for disks 
with strong ordered magnetic field , as fou nd in 2D GRMHD sim- 
ulations by iMcKinnev & Gammid | |2004|) . Of course, in the limit 
that the magnetic field energy density near the black hole exceeds 
the local rest-mass density, a force-free magnetosphere will de- 
velop and deviations from the NT model will become extreme. 
We argued that this corresponds to when the specific magnetic 
flux T > 1 near the disk midplane. Our 1-loop model appears to 
be entering such a phase at late time after accumulation of a sig- 
nificant amount of magnetic flux. Such situations likely produce 
powerful jets that are not observed in black hole x-ray binaries 
in the thermal state. However, transition between the thermal state 
and other states with a stron g power-law component jFender et al.l 
|2004|; iRemillard & McChn tock 2006) may be partially controlled 
by the accumulation of magnetic flux causing the disk midplane (or 
perhaps just the corona) to breach the T ~ 1 barrier. Such a behav- 
ior has been studied in the non-r e lativistic regim e I Naravan et al.l 
I2OO3I : llgumenshchev et alj|2003bl : Ilgumenshchevll2009l) , but more 
work using GRMHD simulations is required to validate the behav- 
ior. 

We also found that the apparently different results obtained 
by NIO were mostly due to measurement and interpretation differ- 
ences. We found that both the disk and the corona-wind-jet con- 
tribute nearly equally to deviations in the total specific angular mo- 
mentum relative to the NT model. However, the corona- wind-jet 
contributes much less to the luminosity than the disk component. 
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Therefore, if one is interested in comparing the luminous portion 
of the disk in the simulations against the NT model, the only fair 
procedure is to consider only the disk gas, i.e., gas within a couple 
of scale heights of the midplane. This is the approach we took in 
this study (also in SOS). NIO on the other hand included the corona- 
wind-jet gas in their calculation of the specific angular momentum. 
The dynamics of the coronal gas differs considerably from the NT 
model. Therefore, while it does not contribute to the luminosity 
of bound gas, it doubles the deviation of the specific angular mo- 
mentum from the NT model. In addition, NIO used a I-loop initial 
field geometry for their work which, as discussed above, further 
enhanced deviations. 
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APPENDIX A: EXAMPLE SOLUTIONS AND SCALINGS 
FOR THE GAMMIE (1999) MODEL 

Table lAll gives representative solutions for the iGammiel j 19991) 
model of a magnetized thin accretion flow. The columns correspond 
to the black hole spin, a; the specific magnetic flux, T; the nominal 
efficiency, e; percent deviation of e from the NT value; the specific 
angular momentum, j; percent deviation of j from NT; and the nor- 
malized rate of change of the dimensionless black hole spin, s (see 
Eg. list. For Y < 0.5 and across all black hole spins, the relative 
change in the specific angular momentum is less than 5% and the 
relative change in the efficiency is less than 9%. For small values 
of T < 1, the deviations of j and e from NT behave systematically 
and one can derive simple fitting functions. For j we find 

iogio[-i)[;]] 

~ 0.79 + 031 (a I M) + 1.601og,oT 
~ (4/5)-Kl/3)(fl/M)-K8/5)log,o 



Table Al. Thin Magnetized Inflow Solutions 



(Al) 
(A2) 



with an L2 error norm of 0.7%, 0.7%, respectively, for the first and 
second relations, while for e we find 



logio [D[e]\ 

1.44 -t- 0.12(fl/M) -t- 1.60 log,o Y 
(3/2) + (l/10)(a/M) + (8/5)log,oY, 



(A3) 
(A4) 



with an L2 error norm of 0.9%, 1%, respectively, for the first and 
second relations. These results indicate that the deviations from NT 
scale as Y*'^ for Y < 1. For Y > 1, the index on Y depends on the 
spin parameter. In the span from Y ~ 0.2 to Y ~ 1, a linear fit across 
all black hole spins gives -£)[;] ~ - 1 -I- 1 1 Y and ~ -4 + 33 Y, 
which are rough, though reasonable looking, fits. 



a 
M 


T 


2 


D\e] 


} 


D[J] 


s 





0.1 


0.0576 


0.709 


3.46 


-0.172 


3.46 





0.2 


0.0584 


2.14 


3.45 


-0.52 


3.45 





0.3 


0.0595 


4.08 


3.43 


-0.991 


3.43 





0.5 


0.0624 


9.17 


3.39 


-2.23 


3.39 





1 


0.0727 


27.1 


3.24 


-6.57 


3.24 





1.5 


0.0859 


50.2 


3.04 


-12.2 


3.04 





6 


0.19 


232 


1.51 


-56.4 


1.51 


0.7 


0.1 


0.105 


1.03 


2.58 


-0.286 


1.33 


0.7 


0.2 


0.107 


3.07 


2.56 


-0.853 


1.31 


0.7 


0.3 


0.11 


5.8 


2.54 


-1.61 


1.3 


0.7 


0.5 


0.117 


12.8 


2.49 


-3.57 


1.26 


0.7 


1 


0.142 


36.7 


2.32 


-10.2 


1.12 


0.7 


1.5 


0.172 


66.3 


2.11 


-18.5 


0.95 


0.7 


6 


0.477 


360 


-0.00714 


-100 


-0.74 


0.9 


0.1 


0.157 


1.17 


2.09 


-0.386 


0.576 


0.9 


0.2 


0.161 


3.37 


2.08 


-1.11 


0.567 


0.9 


0.3 


0.165 


6.29 


2.06 


-2.07 


0.555 


0.9 


0.5 


0.177 


13.7 


2.01 


-4.5 


0.524 


0.9 


1 


0.215 


38.3 


1.84 


-12.6 


0.423 


0.9 


1.5 


0.262 


68.3 


1.63 


-22.5 


0.3 


0.9 


6 


0.845 


443 


-0.958 


-146 


-1.24 


0.98 


0.1 


0.236 


0.949 


1.68 


-0.397 


0.179 


0.98 


0.2 


0.241 


2.86 


1.66 


-1.2 


0.174 


0.98 


0.3 


0.246 


5.36 


1.64 


-2.25 


0.168 


0.98 


0.5 


0.261 


1 1.6 


1.6 


-4.9 


0.152 


0.98 


1 


0.309 


32.2 


1.45 


-13.6 


0.1 


0.98 


1.5 


0.368 


57.1 


1.28 


-24.1 


0.0379 


0.98 


6 


1.21 


416 


-1.27 


-175 


-0.862 


0.998 


0.1 


0.319 


-0.63 


1.4 


0.344 


0.0374 


0.998 


0.2 


0.327 


2.02 


1.38 


-1.11 


0.0342 


0.998 


0.3 


0.332 


3.66 


1.36 


-2 


0.0322 


0.998 


0.5 


0.345 


7.73 


1.33 


-4.22 


0.0273 


0.998 


1 


0.388 


20.9 


1.23 


-11.4 


0.0113 


0.998 


1.5 


0.439 


37 


1.11 


-20.2 


-0.00819 


0.998 


6 


1.19 


272 


-0.675 


-148 


-0.292 



APPENDIX B: INFLOW EQUILIBRIUM TIMESCALE IN 
THE NOVIKOV-THORNE MODEL 

The radius out to which inflow equilibrium is achieved in a given 
time may be estimated by calculating the mean radial velocity v,. 
and then deriving from it a viscous timescale —rjv,-. 

When the flow has achieved steady state, the accretion rate. 



M : 



(Bl) 



is a constant independent of time and position. Here we derive an 
expression for u, corresponding to the general relativistic NT thin 
disk model. In what follows, capital script letters denote standard 
functio ns of r and a (c.f. eqns. (14) and (35) in IPage & Thomd 
( Il974l) ) which appear as relativistic corrections in otherwise New- 
tonian expressions. They reduce to unity in the limit rjM — > oo. 

The vertically-integrated surface density may be defined as 
E = 2/ip, where h is the disk scale-height and p is the rest-mass 
density at the midplane. In equilibrium, density is related to pres- 
sure by 



dp 

— = p X ("acceleration of gravity") 
dz 



■ P- 



Mz &D& 



(B2) 
(B3) 
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the vertically-integrated solution of which is 

h = (p/py'^/\Q.\JiS-'C"-D-^''-8-^'\ 



(B4) 



The pressure may be parameterized in terms of the viscous stress, 
Iff J = Q'P^ which is a known function of r and a: 



M C"-Q 
W = Ihtfi = —Q. . 



The surface density is then 
1 M 



E = 



-Ji^S-^C^'^D--8-'Q. 



Inali^lQ.]^ 

Substituting this in Eq. l IBlt . the radial velocity is 

V, = -a\h/r\'\Q.\r^-^S^C-^'^D^'-6Q-'. 



(B5) 



(B6) 



(B7) 



This result is independent of the exact form of the pressure and 
opacity and so is valid in all regions of the disk. The inflow equi- 
librium time may be estimated as t,^ Ir/v,- 
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